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TABLE OF SYMBOLS 



TEXT E1DV2 


DEFINITION 


A A 


Speed of Sound 


*A 


Denotes the value of a variable at the 
node to the left of a discontinuity, * 
can be any variable name 


Ai/Aj aR 


The ratio of sound speed across a shock. 
A/B (shock moving right) , B/A (shock 
moving left) 




Denotes the value of a variable at the 
node to the right of a discontinuity 


BDRY 


3 denotes left boundary, 2 the right 
boundary 


*BD 


Variable at phantom node 


CNTACT 


3 digit variable denoting contact sur- 
face location, direction of travel, and 
if it crosses a node 


COUNT 


Counter for graphics routines 


CSDIR 


Contact surface direction, 2 to the 
right, 3 to the left 


CSRMN 


Riemann variable change across a contact 
surface 


D2 


Density ratio at first node inside 
boundary 


DARRAY 


Array of density for plotting 


DELAH 


Change in A from I to I+l 


DELAL 


Change in A from I-l to I 


DELQQH 


Change in QQ from I to I+l 


DELQQL 


Change in QQ from I-l to I 


DELRRH 


Change in RR from I to I+l 
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DELRRL 


Change in RR from I-l to I 




DELSH 


Change in S from I to I+l 




DELSL 


Change in S from I-l to I 


6t 


DELT 


Time step 




DELTEX 


The excess time in a time step when the 
shock is exactly at the solid wall 


9t^l 


DELTWL 


The time for the shock to reach the wall 




DELQH 


Change in Q from I to I+l 




DELQL 


Change in Q from I-l to I 


As 


DELX 


Interpolation distance (LMD+DELT) 


P 


DENS 


Density 




DLCD 


Density to the left of the contact 
discontinuity in the exact solution 




DLSH 


Density to the left of the shock in the 
exact solution 




DLTA** 


Prefix which indicates the spatial 
change in ** for one time step 




DQ 


The jump in velocity across the shock 
divided by the sound speed at B (right) 
or A (left) 




DR 


The ratio of the density across a shock, 
B/A (right) , B/A (left) 




DRI 


Initial density ratio across the diaphragm 




EE 


Desired precision for characteristic 
calculations 




E(K) 


Actual error in characteristic slope 
calculation 




EREIMN 


The jump in QQ across the shock calcu- 
lated analytically as a function of w 


Y 


G 


Gamma (ratio of specific heats) 
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GRAPHS 


For graphical output, 0 = none (tabular) , 
1 = plots all variables, 2 = compares 
density with exact solution 


G1 


1/(G-1) 


G2 


2/(G-l) 


H 


1/(N-1) 


HALT 


Terminates program if 1, set by condi- 
tions not coded 


INTEG(K) 


Result of integrating Z (K) 


**INT 


Value of ** interpolated between nodes 
on the current time level 


12 


Number of the node to the right of a 
discontinuity 


JSTOP 


Number of time levels to be calculated 


LBDDR 


Left boundary density ratio 


LBDDRI 


Left boundary density ratio at time zero 


LBDPR 


Left boundary pressure ratio 


LBDPRI 


Left boundary pressure ratio at time zero 


LBDPRS 


Value of 0 denotes constant pressure at 
left boundary, 1 denotes adjustable pres- 
sure at the left boundary 


LBDTR 


Left boundary temperature ratio 


LBDTRI 


Left boundary temperature ratio at time 
zero 


LBNDRY 


Denotes left boundary condition, open or 
closed 


A LMD(K) 


The characteristic trajectories in the 
space-time plane (q+A, q-A , q) 


LNODE 


Array of left most node to be corrected 
in CORRCT 


LWPRES 


Denotes which side of diaphragm has low 
pressure 

1 
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LXX 


Node defining the left interval 




MREIMN 


The measured jump in QQ across the 
shock, from A to B 




N 


Number of Spacial Nodes (odd number) 




ND 


Double precision value of N 




NEW** (I) 


Stored values of ** for the next time 
level 




PARRAY 


Array of pressures for plotting 


Pi/Pj 


PR 


The ratio of the pressure across a shock, 
A/B (right) , B/A (left) 


p 


PRESS 


Pressure 




PRI 


Initial pressure ratio across the diaphragm 




PLTCNT 


Counter for graphics routines 


3* 

3s 


*PRIM(K) 


Suffix which indicates the spatial 
derivative of * at the current time level 




P2 


Pressure ratio at first node inside 
boundary 


q 


Q 


Absolute fluid velocity 




QARRAY 


Array of velocities for plotting 




QLBD 


Initial velocity at left boundary 




QLI 


Initial velocity left of the diaphragm 




QRBD 


Initial velocity at right boundary 




QRI 


Initial velocity right of the diaphragm 


Q 


QQ 


Q+A*S (extended Riemann variable) 




RBDDR 


Right boundary density ratio 




RBDDRI 


Initial right boundary density ratio 




RBDPR 


Right boundary pressure ratio 




RBDPRI 


Initial right boundary pressure ratio 
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RBDPRS 


Value of 0 denotes a constant pressure at 
right boundary, while 1 is for adjustable 
pressure 


RBDTR 


Right boundary temperature ratio 


RBDTRI 


Initial right boundary temperature ratio 


RBNDRY 


Denotes right boundary open or closed 


RNODE 


Array for right most node to be corrected 
in CORRCT 


R RR 


Q-A*S (extended Riemann variable) 


RXX 


Node defining the right interval 


S S 


Entropy 


SAP 


Entropy to the left of the shock for 
flows right or entropy to the right of 
the C.S. for flows left 


SARRAY 


Array of entropy for plotting 


SAVG 


Average entropy 


SBP 


Entropy to the right of the C.S. for 
flows right or entropy to the left of 
the shock for flows left 


SHKDIR 


Shock direction of travel, 3 to the left, 
2 to the right 


SHOCK 


3 digit variable denoting shock location, 
direction of travel, and if it crosses a 
node 


SIGMA 


Spatial location of discontinuities 
Sigma (L,J) where L indicates the type of 
discontinuity and J indicates the time 
level: 1 — current level, 2 — level being 
calculated 


SK 


Integer that denotes relative location of 
shock near boundaries 


SKIP 


Variable which indicates how many time 
steps between calls to output routines 


**STEP 


The change in time of ** at a node used 
to step up to the next time level 
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t 


T 


Time since initial conditions 


T 


TEMP 


Temperature 




TRI 


Initial temp, ratio across the diaphragm 




TS 


Time for shock to travel one interval 




T2 


Temperature ratio at first node inside 
boundary 




UA 


Velocity relative to the shock, left side 




UB 


Velocity relative to the shock, right 
side 




VHEAD 


Velocity of the head of the expansion 
wave for the exact solution 




VCDE 


Velocity of the contact discontinuity 
for the exact solution 


Vs 


VS 


The shock speed (positive right, 
negative left) 




VSE 


Velocity of the shock for the exact 
solution 




VTAIL 


Velocity of the tail of the expansion' 
wave for the exact solution 


w 


W 


Mach no. relative to a standard shock 


w 




Vector of the principle variables, Q,R,S 


X 


X 


Location in spacial plane (I-1)*H 




XARRAY 


Array of spatial positions for plotting 




XEXACT 


Array of six X values for the exact 
solution 




XINIT 


Initial position of discontinuity for 
exact solution plotting 




X2 


Location of node to right of discont. 
along the spacial axis 




Y 


(N+D/2 




YEXACT 


Array of six density values for the 



exact solution 
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7 . 



Z(K) 



A 



Vector of the right hand sides of the 
governing equations 

Small spatial change 



6 



Small change with respect to time 



e,4> 



Flow angles with respect to reference 
coordinate planes 
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I. 



INTRODUCTION 



The investigation of wave rotors^ and their possible 
application in gas turbine engines at the Naval Postgraduate 
School's (NPS) Turbopropulsion Laboratory prompted the 
development of a one-dimensional unsteady flow computer 
code. The QAZID method developed by Verhoff [Ref. 1] was 
chosen over other methods [Ref. 2] for reasons which 
included the following: 

1. The method is based on the use of characteristics. 
Such methods can model wave propagation accurately. 

2. The use of a natural streamline coordinate system 
eases the difficult task of computing with two and 
three dimensional grids. 

3. The equations are written in a form which allows a 
straightforward extension to viscous flows. 

Salacka [Ref. 2] verified the development of the equations 

and implemented a one-dimensional Euler solution of the 



^The wave rotor consists of simple tube-like passages 
along and around the periphery of a drum which rotates 
between end walls containing inlet and outlet (partial 
admission) gas ports. Compression and expansion processes 
are arranged to occur in a cyclically periodic unsteady 
process within the rotor such that a high pressure driver 
gas is used to compress a low pressure driven gas. The 
compressed driven gas is led from the outlet port to the 
combustor and brought back into the rotor as the driver. 
The expanded driver gas exits an outlet port. Thus the 
rotor functions as both turbine and compressor with direct 
gas-gas energy exchange through unsteady wave processes. 
The technology of such devices requires the design of 
appropriate cycles using one-dimensional calculations, and 
analysis of the multi-dimensional unsteady flows such as 
occur during port opening and closing. 
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shock tube problem in the EULERl code. The EULERl code was 
limited to high pressure set on the left, with flow and 
shock velocities to the right. The shock wave was tracked 
and corrected for in the flow, but the contact discontinuity 
was not, and the expansion wave was not tracked. The code 
also did not treat end boundary conditions. 

The present work extended the EULERl code of Salacka to 
become the "Euler ID Version 2" or E1DV2 code which: 

1) can handle both right and left traveling flow and 
discontinuities 

2) implements tracking and correction of the contact 
discontinuity 

3) implements tracking of the expansion wave 

4) models closed wall and open end boundary conditions 

5) models shock and contact surfaces within a grid 
interval . 

Section II and Appendix A describe the analysis underlying 
the revisions. The Fortran code, E1DV2, is detailed in 

Section III and Appendix C and graphical results of four 
test cases to demonstrate the new features of the code are 
given in Section IV. A discussion of these results and 

limitations of the present code follow in Section V. 

Conclusions and recommendations are made in Section VI. 

Appendix B contains instructions to operate the code on the 
NFS computer and the FORTRAN listing is included as Appendix 
D. 
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II. ANALYTICAL AND COMPUTATIONAL APPROACH 



A. QAZID DESCRIPTION 

The QAZID analysis and numerical solution scheme is an 
explicit, non-conservative method that uses extended Riemann 
variables to model the multi-dimensional Euler equations in 
a natural streamline coordinate system. Table I lists a 
summary of the applicable equations which are derived in 
detail in [Ref. 2]. 

1 . The Coordinate System 

Figure 2.1 shows the natural streamline coordinate 
system (s,n,m) relative to a fixed rectangular cartesian 
system (x,y,z). The right-hand orthogonal system moves in 
curvilinear translation along a streamline. The "s" 
direction is measured in the direction of the flow, tangent 
to the streamline. The "n" direction is perpendicular to 
the s direction in a plane perpendicular to the x-z plane. 
The "m" direction is normal to the plane containing the n 
and s directions. The angle is the angle between the s-n 
and x-y planes, and the angle 9 is the angle between the 
velocity vector and the x-z plane in the s-n plane. In the 
one-dimensional problem treated here, the s direction is 
such that velocity to the right is positive, and to the left 
is negative. [Ref. l:p. 2] 
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TABLE I 



SUMMARY OF THE EULER EQUATIONS 

Definitions 



Speed of sound 




A = /fP/p 


(I.l) 


Modified entropy 








change 




dS - - ^ 

yT 


(1.2) 


so that 










S 


■ ^ ■ -(7pr)C2Y - ln(P/p^)] 


(1.3) 


Extended Riemann 








variables 




Q = q + As 


(1.4) 






R = q - As 


(1.5) 



Euler Equations 



3t 



+ (q+A)|| 



- - ;^) [|j(q - ;^)J 

- ^ 9 s' 



( 1 . 6 ) 



3R ^ , ,.3R 



ds 



+ :i^qAs(||+oos e|i] 



3t 



3S 
^ 3s 



3t 



rr M 

^ 3s 



3t 



+ q 



3s 



^ 3 In P 

Yq 8n 

t? 3 In P 
yq OOS0 3m 



(1.7) 

( 1 . 8 ) 
,(I-9) 
(I. 10) 
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Figure 2 . 1 The Natural Coordinate System 



2 . Extended Riemann Variables 

The "Riemann Variables" in most publications are 
defined as 



and (2.1) 

R2 = <J + (tIt) ^ 

where q is the velocity magnitude, A is the speed of sound, 
and y is the ratio of the specific heat at constant pressure 
to the specific heat at constant volume for a particular gas 
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[Ref. 3:p. 5], Verhoff modified the definition to obtain 

"Extended Riemann Variables," defined as 

Q = q + AS (1.4) 

R = q - AS (1.5) 

where S is the modified entropy given in Table I [Ref. l:p. 

2 ]. 

3 . The Governing Equations 

The Euler equations in Table I were developed in 
detail [Ref. 2: Appendix A] from the equations for 

conservation of mass, momentum, and energy for a control 
volume. The underlying assumptions comprise inviscid flow, 
negligible body forces, no heat transfer, and a perfect gas. 
These assumptions and the definition of modified entropy 
were used in transforming the equations into a form suitable 
for solution using the method of characteristics. 

B. CHARACTERISTIC SOLUTION 

Equations (1.6), (1.7), (1.8), (1.9) and (I. 10) are a 

system of slightly coupled quasi-one-dimensional partial 
differential equations (PDE) with eigenvalues q+A , q, and q- 
A. These equations can be rewritten along the characteris- 
tics in the time-space domain to become ordinary 
differential equations (ODE) which are solved by quadrature 
and interpolation. The propagation of each characteristic 
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curve along its particular trajectory is expressed by the 
respective ODE. [Ref. l:p. 1] 

Each in the system of equations is in the form 



9w 9w 
9t ^ 9s 



z 



( 2 . 2 ) 



If w = w(s,t) it is shown by Salacka [Ref. 2:pp. 20-23] that 



A 



ds 

dt 



(2.3) 



and 



where A describes the characteristic curve along which w 
changes . 

In this section the solution of the 3D Euler equations 
given in Table I is described. The computer code developed 
in the present work, however, was for the ID case wherein 
equations (1.9) and (1. 10) are not required, and only the 
three remaining equations (1.6), (1.7), and (1.8) are 

considered. 

Equations (1.6) to (1.8) for one-dimensional flow may be 
expressed in matrix form by setting 
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- 

Q 




q+A 


O 

O 

1 


w = 


R 


[X] = 


0 


0 

< 

1 

cr 




S 




0 


0 q 













and 



2 = 



Y_i) [q 



-^) [q +;^]s 



Then equation (2.2) becomes, for the equation set. 



i " i 



z . 



Figure 2.2 shows the characteristic curve in the space- 
time domain between two nodes within the computational mesh. 
Equation (2.2) has a solution 

_ _ t4dt _ 

(Sw = - Aw + / z dt (2.5) 
t 

where 6w is the change due to time at a fixed location, and 
Aw is the change due to displacement at a fixed location, 
along the characteristic curve. The value of Aw is 
calculated by interpolation through an iterative scheme 
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5w = w_ - w 
B s . 



Aw = w - w, 

S. A 

Figure 2.2 Characteristic Curve Within Computational 
Mesh 



using initial guesses for the characteristic slope, X , from 
values at time level t. The line integral is approximated 
by 



t+dt 

/ 2 dt 




= 2(As/X) 



S. 

1 — 



S.-AS 
1 



/ 



f(S.-(S,-AS)) 

AX X 



= z(5t) 



( 2 . 6 ) 
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The spatial derivatives of q and A in "z" are estimated from 
values at A and s^. Finally, 3w is calculated, and each 
node is updated to the next time level. 

The Courant-Friedrichs-Lewy (CFL) convergence condition 
requires the domain of dependence of the numerical approxi- 
mation to include the domain of dependence of the differen- 
tial equation [Ref. 4:pp. 336-337]. Figure 2.3 shows that 




Figure 2-3 Computational Grid for the QAZID Method 

the computational grid for the QAZID method satisfies this 
condition by having interpolated initial data points (points 
1, 2, and 3) in between previous solution nodes (points 4, 
5, 6) . The q+A characteristic is estimated from values at 

point 4, q characteristic from values at point 5, and q-A 
characteristic from values at point 6. The use of a maximum 
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time step keeps the domain of dependence of the numerical 
scheme just outside that of the actual equations, aiding 
convergence. Unlike finite difference schemes that are 
conservative and require a numerical diffusion or viscosity 
to handle oscillations or smearing, this method is stable 
without numerical smoothing or artificial dissipation [Ref. 
5:pp. 562-583]. However, the equations do not account for 
the discontinuities across shocks and contact surfaces 
properly. Thus after each time step a one-point correction 
technique is used to correct for shocks and contact 
surfaces. 

C. DISCONTINUITIES 

1 . Exoansion/Rarefaction Waves 

The head and tail of the expansion wave travel along 
the characteristic as gradient discontinuities. Flow 
velocity, speed of sound, pressure, and density are continu- 
ous across a gradient discontinuity but the spatial deriva- 
tives are discontinuous. Entropy remains constant through 
the rarefaction wave. The head of an expansion wave moves 
into undisturbed flow with a speed q+A, while the tail of 
the expansion wave has a velocity q-A. No special treat- 
ment of these waves is required; the characteristics method 
accurately tracks and models their influence. Expansion 
waves are generated when two shocks collide, a diaphragm is 
burst, a contact surface and shock intersect, a shock leaves 
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an open boundary, or an open boundary has a pressure lower 
than the pressure inside the tube. [Ref. 3:pp, 13-15] 

2 . Shocks 

Shocks are created when a diaphragm is burst, 
compression waves generated at an open boundary coalesce 
into a wavefront with sufficient strength, or when a shock 
and contact surface collide. Pressure, velocity, density, 
and entropy are discontinuous across a shock [Ref. 6:pp. 30- 
31]. The equations of motion (I. 6) -(I. 10) are not correct 
for calculating changes through shocks. A method to correct 
for a shock is to use the ratio of the jump in the extended 
Riemann variable across the shock to the speed of sound 
downstream of the shock to find the incoming Mach Number 
relative to the shock (w) [Ref. l;p. 3]. Figure 2.4 shows 
the case of a shock headed to the right. The shock velocity 
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Figure 2.4 Shock Wave with High Pressure on Left 
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(Vg) is imposed on the flow to produce a stationary shock 
condition. Flow variables are then corrected using normal 
shock relations. Figure 2.4 applies to a condition where 
high pressure is on the left. The shock is moving to the 
right at Vg into flow with velocity qg. Table II lists the 
equations developed by Salacka [Ref. 2] for relating the 
extended Riemann variable change to Mach Number. Velocities 
are non-dimens ionali zed by the speed of sound downstream of 
the shock direction, Ag, and pressures and densities by the 
respective values downstream. 

Figure 2.5 depicts equation (II. 2) over a range of 
Mach Numbers from 1.0 to 4.0. The curve is approximated by 
the quadratic equation 






'B 



-2.7574 + 3.1573W - 0.2863w2 



(2.7) 



If the high pressure is to the right the shock moves 
as Fig. 2.6 illustrates, and the governing equations are as 
listed in Table III. Again, downstream conditions are used 
to non-dimensionalize velocity, pressure, density, and the 
Riemann variable change. The development of the equations 
in Table III is given in Appendix A. A plot of equation 
(III. 2) over a range of w from 1.0 to 4.0 is shown in Fig. 
2.7. The curve is approximated by the quadratic equation 



v-\ 

~xr 



2.7574 - 3.15732W + 0.2863w2 



( 2 . 8 ) 
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TABLE II 



SHOCK EQUATIONS FOR HIGH PRESSURE ON LEFT 



Definitions 



Relative Incoming 
Mach Number: w 






(qB-Vs) 

"b 



(II. 1) 



Equations 

Q Riemann 
Variable Change: 



^a”^b 



2(w^-l) 

(Y+1)w 



2 . /a 






Ar 1 



n rr 2Y 2 
■] In [ (—-r W 



'Y+1 



- (^)) ■^2) ^1 
(Y+l)w 



Speed of Sound Ratio: 



(II. 2) 



% = 7- j; ^ .- [2(Y-l) [1 + ^ w^] [-^ - 1]]^"^^ 

(Y+1)w " ' 2 Y~1 



(II. 3) 



Pressure Ratio: 



^ = JL _ III 

Pg Y+1 Y+1 

Density Ratio: 

^ = (Y+l)w^ 

•^B (y-1)w^+2 



(II. 4) 



(II. 5) 
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TABLE II (CONTINUED) 



Velocity Change; 






2(w^-l) 

(Y+1)w 



Entropy Change: 



V^B = 



-( L_)in[2DlZl±l] 



, 1 

^ (Y+1)w^ 



(II. 6) 



(II. 7) 
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Figure 2.5 Q, Extended Riemann Variable Change with Mach Number 
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Figure 2.6 Shock Wave with High Pressure on Right 
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TABLE III 



SHOCK EQUATIONS FOR HIGH PRESSURE ON RIGHT 



Definitions 



Relative Incoming 
Mach Number: 



w 





(III-l) 



Equations 

R Riemann 
Variable Change: 






A 



2 (1-w^) 
(Y+l)w 



2 ^ 
A 



B r 1 m r / 2y , 2 

'•^Y+1 



- (ari))( (Y-i)« -^2 )Yj 



Y+1 



(y+l)w 



(III. 2) 



Speed of Sound Ratio: 



r = V [2 (y-l) [1 + W^] [-^ - 1 ] ] 

K (y+l)w 2 '■y-l 



1/2 



(III. 3) 



Pressure Ratio: 



^ = Jl. 2 _ y^i 
y +1 y+i 



Density Ratio: 



% = (y+i)w^ 

^B (y-l)w^ +2 



(III. 4) 



(III. 5) 
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TABLE III (CONTINUED) 



Velocity Change: 



^ 2(l-w^) 

(y+1)w 



(III. 6) 



Entropy Change: 



o 
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Figure 2.7 Extended Riemann Variable Change with Mach Number, High 
Pressure on Left- 



An examination of equations (II. 2) and (III. 2) 



reveals that 







This is expected since in flow to the right, the Q extended 
Riemann variable is associated with the q+A characteristic. 
Also, since velocity is defined as positive to the right, q 
is positive. But for flow to the left, velocity is defined 
as negative, thus q is negative. Moretti [Ref. 3] showed 
that the downstream running Riemann variable 
(q+A characteristic) has a lower magnitude of discontinuity 
across a normal shock than the upstream running Riemann 
variable (q-A characteris-tic) , and is thus preferable. In 
flow to the left, the Riemann variable associated with 
downstream is the R Riemann variable, though the character- 
istic associated with the R Riemann variable is q-A , in flow 
to the left q is negative and thus 



Therefore, in fluid traveling left the R Riemann variable is 
used to find Mach Number, and in fluid traveling right the Q 
Riemann variable is used. 



q-A = - ( 1 q 1 + A ) 



(2.9) 
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If the shock is between two nodes with no contact 



surface within the same interval, the method to calculate 
the Mach Number, w, is as follows: 

1) The change in the appropriate extended Riemann varia- 
ble across the interval is measured; R variable for 
flow left, and Q variable for flow right. Call this 
change , Am. 

2) Use Am in equation (2.7) or equation (2.3) for the 
appropriate flow to calculate w. 

3) With equation (II. 2) or equation (II. 3), for flow 

right or left respectively, use w to calculate the 

analytical Riemann variable change, Ae. 

4) If this exact value of Ae is not equal to Am then 

calculate a new guess, Ag from 

Agi+i = Agj_ + (Am - Ae) (2.10) 

and then repeat steps 2 to 4 until Ae equals Am. 

Thus the correct value of w is determined, and shock 
corrections from the normal shock relations are valid. Vg 
is determined from equation (II. 1) or equation (III.l) 

depending on flow direction. Flow direction will be 
initially determined by the side with high pressure. High 
pressure on the left means flow will travel to the right. 
Likewise, high pressure on the right means flow will travel 
to the left. 

3 . Contact Surface 

A contact surface is a boundary between regions of 
flow which are different in composition or which have 
undergone different thermodynamic histories. Thus a contact 
surface is formed when the diaphragm is burst in a shock 
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tube separating the gas initially in the driver from the gas 
initially in the driven tube [Ref. 6:pp. 23-29]. In a wave 
rotor, a contact surface would separate the combustor gases 
from the cooler inlet air. Contact surfaces are also caused 
by two shocks of opposite families colliding or a shock over 
taking a shock of the same family [Ref. 3:pp. 15-18], 

Figure 2.8 illustrates the changes in physical 
properties through the contact surface in the shock tube 
problem. The gas to the right has been compressed and 
heated by the shock wave, with that to the left cooled and 
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Figure 2.8 Behavior ^of Physical Properties Through 
a Contact Surface 



40 



expanded. The density, speed of sound, and modified entropy 
are discontinuous across the contact surface. 

Since pressure is constant across the contact 
surface, using the perfect gas equation of state and the 
first law of thermodynamics it is shown in Appendix A for 
the shock tube problem that 



(y-1) 



V"b 



= exp 






( 2 . 11 ) 



Since the flow behind the contact surface to the 
tail of the rarefaction wave is isentropic, the value of Sg 
will be known. The value of is the value of entropy 
behind the shock. Subtracting equations (1.4) and (1.5) 
from each other and dividing by entropy Sg, Ag can be 
determined. Thus the unknown speed of sound A a obtained 
from equation (2.11). With qg equal to and Sp^ the 
Riemann variables just behind the contact surface are calcu- 
lated from equations (1.4) and (1.5). For a contact surface 
traveling to the left it is only necessary to interchange 
the subscripts. The velocity of the contact surface is 
easily determined since it moves at velocity q along the q 
characteristic line associated with S. [Ref. 3:p. 16] 

4 . Contact Surface/Shock Interaction 

When the shock and contact surface are within the 
same interval the calculation of the shock incoming Mach 
Number, w, must be modified. Figure 2.9 shows a plot of the 
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Entropy 




Figure 2.9 Modified Entropy Distribution for Shock 
and Contact Surface Traveling Right 



modified entropy change across the interval. The use of the 
Riemann values at the nodes to estimate the Riemann variable 
change across the shock would be incorrect. The correct 
Riemann variable change would be, for flow to the right. 



ASK = (2.12) 



The correct Mach Number, w, is determined using the 
technique of Moretti modified for QAZID [Ref. 3:pp. 23-24]. 

Referring to Fig. 2.9, the technique for flow to the 
right is as follows: 

1) Calculate w as if no contact surface were in the 
interval. With values at node i +1 known, determine 
S 3 . Then use w and S^ in equation (II. 7) to solve for 
S;^. This is the initial estimate of S 3 . Let S 3 ' = S 3 . 

2) The Q Riemann Variable change across a contact surface 
is shown in Appendix A to be given by 
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Acs 



(2.13) 




3) Define 



Am 




(2.14) 



then 



Am = ASK + ACS 



(2.15) 



can be solved to obtain ASK. 

4) Using ASK as Am in the shock iteration scheme given in 
Section II. C. 2 calculate a new w. 

5) With this new w and S 3 use equation (II. 7) to obtain a 
new S;^. 

6 ) Compare Sji^ with S 3 . If they are not equal to within 
an acceptable error, calculate a new estimate for Sc 
using 



Iterate until convergence is achieved. 

This will result in the proper w for a condition 
with shock and contact surface interacting. 

For flow to the left, it is only necessary to inter- 
change the subscripts A and B, and use 




(2.17) 
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for the extended Riemann Variable change. Figure 2.10 
illustrates the left traveling condition. 




Figure 2.10 Modified Entropy Distribution for Shock 
and Contact Surface Traveling Left 

D. BOUNDARY CONDITIONS 
1. Open Boundary 

At the open boundary, a reference pressure ratio, 
Pj.gf = Poo/Pa specified, where P„ is the pressure to be 
held at the boundary, and Pj^ is the pressure just inside the 
tube. Only the case where P^ < ?a considered. This 

means that an outflow condition at the boundary will result 
when < P;^, with an expansion wave traveling in from the 
boundary. Thus locally at the boundary, conditions are 
isentropic. 

The computational grid for the left boundary is 
shown in Fig. 2.11. A phantom node, L, is used outside the. 
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left grid boundary 
B 




Figure 2.11 Left Boundary Computational Grid 

computational mesh to enable simple enforcement of constant 
pressure and entropy at the mesh boundary node, 1. 

The approach is to determine the required variables 
at node L, then transfer these values onto node 1, and vary 
q at the boundary to meet the required boundary conditions. 
First, using P^ef node L and S at node 2 in equation 
(1.3) rewritten in the form 

(s--^)(-y(y-D) -1/y 

Pr = {.(1/Pj.ef) (exp )) (2.18) 
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Pl^, the density ratio at the boundary, is calculated. Using 
the perfect gas equation of state in equation (2.18) the 
temperature ratio, Tj^ is given by 






(Y-1) 

(exp ^ ) 



(2.19) 



and can be calculated once is known. With an initial 

velocity at the boundary, qj^, assumed the Riemann variables 
R and Q are determined for node L using 



R = 




( 2 . 20 ) 



and 



Q =. <3b + S (2.21) 

Subtracting equation (2.20) from (2.21), and 
dividing by twice S yields A at node L, where A = /y RqT. 
The values of R, Q, A, qt), and S are now known at node L. 
These values are assigned to node 1, and the QAZID algorithm 
is applied to the boundary node as if it were an interior 
node. The result is an estimate of the new values for R, Q, 
and S at node 1. These are used to obtain a new estimate of 
the pressure ratio at node 1, call it Prn* ^RN compared 
with Pref' value of S with the old value of S at 
node 1. A new qj^ is calculated by solving equations (2.20) 
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and (2.21) simultaneously. The process is repeated until 
entropy and pressure remain constant at the boundary. For 
the right-hand boundary, Fig. 2.12 shows the computational 
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Figure 2.12 Right Boundary Computational Grid 

grid with the phantom node, r, to the right of the mesh 
boundary node, n. The procedure for the right boundary is 
the same as for the left boundary. 

When a shock travels across the open boundary, the 
boundary is first corrected using the method described in 
Section II. C. 2. Subsequently, the pressure is returned to 
give the constant value specified for P^ef' using the above 
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procedure. The result will be an expansion wave traveling 
inward . 

For supersonic flow, all of the variables are deter- 
mined from values at the interior nodes [Ref. l:p. 3]. 

2 . Closed Boundary 

The solid boundary is imposed by setting the 
velocity at the boundary, q, equal to zero. The same 
computational grid is used as that for an open boundary, but 
a different technique is used in assigning values to the 
phantom node. Referring to Fig. 2.11, the velocity at node 
2, q 2 , is known. By setting 

qL = - q2 (2.22) 

the wave impacting the boundary will be met by a wave of 
equal velocity but opposite in direction and will result in 
node 1 appearing as a solid boundary. Further, to 

facilitate the QAZID method, set 



Rl = 


- Q2 


(2.23) 


Ql = 


IR 2 I 


(2.24) 


Sl = 


^2 


(2.25) 


'^L = 


• 2'A 


(2.26) 
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This will enable the computation of the new boundary node 1 
values as an interior point with 

qi = 0 

The right boundary is handled in the same way. 

Figure 2.13 illustrates the sequence of events when 
a shock reflects off a solid boundary on the right [Ref. 7: 
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Figure 2.13 Shock Reflecting at Solid Boundary 

pp. 172-195]. Since the velocity behind the reflected 
shock, qg, is zero, and the velocity in front of the shock, 
q;^, is known, the velocity gradient, Aq, over the shock is, 
in non-dimens ionalized form, given by 
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Aq 



(<3B “ qA)/^A 



(2.27) 



Rearranging equation (III. 6 ), and using Y = 1.4 it 
can be shown that 

w = \/l + 0.36(Aq)2 - 0.6Aq (2.28) 

This is an exact analytical value for w based on stationary 
normal shock relations. The speed of sound ratio, ^b/'^A' 
calculated from equation (III. 3) and 

B “ (2.29) 

Entropy behind the shock, S 3 , is determined from equation 
(III. 7) since is known. Then the extended Riemann 

variables, R and Q, behind the reflected shock are 

calculated from equations (1.4) and (1.5) respectively. 

For a shock reflected from the left boundary, the 

subscripts A and B are interchanged and using equation 
(II. 6 ) the Mach Number is given by 

w =yi + 0.36(Aq)2 + 0.6Aq (2.30) 

The equations from Table II are used instead of 

those from Table III (used in the case of reflection from 
the right boundary) , since the shock is now traveling from 
left to right. 
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III. FORTRAN PROGRAM E1DV2 



A. GENERAL DESCRIPTION 

Program E1DV2 is a Fortran computer program which 
calculates 1-dimensional unsteady flow based on the QAZID 
formulation and method of solution of the Euler equations of 
motion. The program has provisions for tracking in time the 
location of shock waves, contact surfaces and head and tail 
of the expansion waves. Its intended application to wave- 
rotor design and analysis explains the present definition of 
node points and boundary conditions. In its current form, 
the code describes flow in a tube initiated by a diaphragm 
bursting. The location of the diaphragm, whether the flow 
is to the left or the right, and whether the ends of the 
tube are closed or open, can be varied. 

The E1DV2 program is written in FORTRAN 77, and was 
developed using the IBM 3 03 3 System 370 computer. The use 
of DISSPLA subroutines in plotting results requires 
compiling and executing using VS FORTRAN procedures. 
Appendix B describes the procedures for running the E1DV2 
program on the Naval Postgraduate School IBM computer 
system. A listing of the code, which contains its own 
description in comment statements, is given in Appendix D. 
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The code incorporates: 

1) first order accuracy 

2) doxible precision numerics, except for single precision 
in graphical output 

3) linear interpolation and extrapolation algorithms 

4) quadratic polynomial approximations for curve fitting 

5) non-dimensional variable input and output 

6) structured subroutine format. 

Table IV lists all the subroutines called from the main 
program and brief descriptions of their purposes. The 
extensive use of subroutines was intended to allow modifica- 
tions and extensions of the code without major 

restructuring . 

B. THE MAIN PROGRAM 

The main program flowchart is illustrated in Appendix C, 
Fig. C.l. The conventions and a comprehensive list of 
variables are listed in the beginning of the program. The 
number of grid nodes is set by the user, and must be an odd 
number. Arrays are dimensioned to the number of nodes in 
the main program prior to compiling and executing. 

Initial values for temperature, pressure, and density 
ratios at the diaphragm and boundaries are entered by 
editing. Also, the initial velocity distribution on each 
side of the diaphragm is set, and as well as the ratio of 
specific heats. The side with the high pressure is identi- 
fied. Boundaries are set either closed or open. In the 
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TABLE IV 



SUBROUTINES 



TIME 

TRAK 

SWEEP 

CONDI 

COND2 



COND3 



COND4 



COND5 



COND6 



COND7 



COND7N 



Calculates minimum time step 

Tracks discontinuities, calculates shock mach 
number, w 

Advances node values to next time step 

Calculates interpolated values of Riemann 
variables, 3q/8s, and 9A/ 3s when no shock 
or contact surface within an interval 

Calculates interpolated values of Riemann 
variable, 3q/3s, and 3A/3 s when discontinuity 
within interval to right of node or flow is 
supersonic to the right 

Calculates interpolated values of Riemann 
variable, 3q/3s, and 3A/3s when discontinuity 
within interval to left of node or flow is 
supersonic to the left 

Called when node is jumped by discontinuity 
to load node information into array for use 
in subroutine CORRCT 

Calculates interpolated values of Riemann 
variable, 3q/3s, and 3A/3s when discontinuity 
on either side of node and flow is 
supersonic 

Would contain subroutine to advance node 
when contact surface is in front of shock 
after intersecting and traveling in same 
direction. Not currently in use 

Would contain subroutine to solve precisely 
when and where shock and contact surface 
would intersect within an interval . Not 
currently in use 

Would contain subroutine to solve precisely 
when and where shock and contact surface would 
intersect if either were jumping a node simul- 
taneously. Not currently in use 
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S 



TABLE IV (CONTINUED) 



COND7S 


Would contain the subroutine to solve the 
Rieiaann problem of shock and contact surface 
intersecting. Not currently in use 


COND8 


Would contain subroutine to advance node after 
shock and contact surface have crossed and are 
moving apart. Not currently in use 


CORRCT 


Advances nodes jumped by discontinuity to next 
time step 


DELTAX 


Determines relative location of discontinuity 
within an interval 


SKJUMP 


Calculates variable change across a normal 
stationary shock 


CSJUMP 


Calculates variable change across a contact 
surface 


EXTRAP 


Extrapolates entered values 


INTERP 


Interpolates entered values 


DBURST 


Calculates initial Riemann values at node 
where diaphragm is located 


BONDRY 


Calculates values to update left and right 
boundary 


SRFLCT 


Calculates new shock parameters when shock 
reflects off solid boundary 


BBDRY 


Alternate interpolating scheme near boundary 


BORDER 


Sets graphical borders 


PLOT 


Plots q, S, P, and p 


EXACT 


Plots exact p versus calculated p 


LIST 


Lists calculated values in files 8, 9, 10 



open case, either constant pressure, or an adjustable 
pressure is specified. The latter case in effect extends 
the tube, and allows discontinuities to disappear out of the 
tube without creating expansion or compression waves. Exact 
values for velocities of expansion wave, shock, and contact 
surface, and density ratio are specified by the user in the 
input section. 

Output is controlled by setting the variable GRAPHS to 
either 0, 1, or 2 . A zero creates three files which contain 
data calculated at the time the call to subroutine LIST is 
made. A value of 1 causes plots of pressure, density, 
velocity, and entropy distributions to be created. When 
GRAPHS equals 2 a plot of the exact density distribution 
along with the calculated density distribution is made. 
Calls to output routines are determined by the parameter, 
SKIP, the number of time steps between calls, which is 
specified by the user. 

Program termination can occur for several reasons. 
First, the program terminates when a maximum number of time 
steps, JSTOP, is reached. JSTOP is specified by the user. 
Second, the program terminates if, during execution, any of 
the following conditions are met: 

1) The shock intersects with the contact surface 

2) The pressure at any open boundary is higher than the 
pressure at the first node inside the boundary 

3) The shock exits an open boundary. 
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In these cases a message is displayed on the terminal 
describing the reason for termination. 

C. THE SUBROUTINE PROGRAMS 

1. The **DBURST” Routine 

Figure C.2 in Appendix C shows the flowchart for 
"DBURST" . The assumption that a shock and contact surface 
are formed immediately causes oscillating numerical solu- 
tions that dampen within five to ten time steps. By 
mathematically "bursting” the diaphragm prior to time zero, 
a solution at the node where the diaphragm is located can be 
determined at time zero. The result is a more accurate 
representation of the wave structure and a dampening of a 
reduced transient solution within three time steps. 
"DBURST" calculates the Riemann variables that would exist 
behind a shock and contact surface that have moved an 
infinitesimal amount after the diaphragm burst, and assigns 
them to the node where the diaphragm was situated. 

2. The "TIME" Routine 

The "TIME" subroutine was not changed from that 
developed by Salacka [Ref. 2:p. 35]. The purpose is to 

compute the maximum allowable time step but ensuring that 
all characteristic trajectories over the time step are 
within one interval between nodes. The time step is 

determined, from 
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DELT 



H/ABS[Q4A ] 



(3.1) 



at every node, and the minimum time step is used. 

3 . The ”TRAK” Routine 

This routine computes the new locations of the 
shock, contact surface, and head and tail of the expansion 
wave after the time step, DELT, calculated in "TIME." In 
Appendix C, Fig. C.3 shows the flowchart for the "TRAK" 
subroutine. The shock velocity, Vg, incoming Mach Number, 

w, and pressure, density, and sonic velocity ratios PR, DR, 
and AR respectively are determined. The equations for flow 
left and right are the same, except for DQ, the velocity 
gradient, which is different only by sign. Thus by proper 
bookkeeping the same equations are coded once and used in 
either direction. The initial direction of the shock and 
contact surface are set from code that interprets which side 
has the high pressure at time zero. The contact surface 
speed is determined after the shock speed is established. 
The situation at time zero is handled as in "DBURST" to 
ensure that the shock/contact surface interaction is 
properly modeled. Since there is a certain transient 
solution that decays after three time steps, tracking of the 
expansion wave is not initiated until then. 

The head of the expansion wave travels at a velocity 
corresponding to q+A , and the tail with speed q-A. In the 
case of flow traveling to the left, this is reversed. 

♦ 
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4 . The »SWEEP" Routine 

The "SWEEP" subroutine solves equation (2.5) and 
updates the values of variables at the nodes, including the 
boundary nodes. Appendix C, Fig. C.4 shows the "SWEEP" 
flowchart. Only those nodes which have been crossed by a 
discontinuity during the time step are updated by the 
"CORRCT" subroutine. Section II. B describes the theory 
coded into "SWEEP". The routine begins at the left boundary 
node and progresses to the right boundary node, one node at 
a time. Figure 3.1 illustrates the overall algorithm 
applying the QAZID method. At the boundary nodes, the 




Figure 3.1 Overall Algorithm for QAZID Method 
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"BONDRY” subroutine is called and returns the updated values 
at those nodes. In between, the correct and most effective 
algorithm to update the node is determined. These 

algorithms are coded into eight different condition 
subroutines. These routines calculate the interpolated 
Riemann variables associated with the three characteristics 
(q+A, q, and q-A ) , plus the spatial derivatives at 9q/3s and 
9A/3S* This allows "SWEEP" to calculate Aw and z. The 
integral of ~z is then determined, and the update for the 
variables of that node in 3w are computed. 9w is stored 
until the entire mesh has been swept. Then the variable 
arrays (w) are updated to the next time interval. 

To choose the proper condition routine, the 
following information about conditions in the interval on 
either side of the node are expressed in two 3 digit integer 
variables, SHOCK and CNTACT. The value of these variables 
provide a code for the following information: 

SHOCK — The existence of a shock within an interval, and, 
if so, the direction of travel, and if the shock 
will cross the node in front of it within the 
current time step. 

CNTACT-The existence of a contact surface within an 
interval, and, if so, the direction of travel, and 
if the contact surface will cross the node in front 
of it within the current time step. 

Figure 3.2 illustrates the regions around a node. Table V 

lists the values that SHOCK and CNTACT may have, and their 

meaning. Figures 3.3, 3.4, and 3.5 illustrate examples of 

what different SHOCK and CNTACT values mean. At each node 
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2 Interval and Node Description used in "SWEEP” 
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3.3 Schematic of SHOCK = 331 and CNTACT = 232 
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TABLE V 



VALUES OF PARAMETERS SHOCK AND CNTACT 



SHOCK 



Located 

in interval on 


Direction 
of travel 


Crosses node in 
oath within time ; 


No Shock 


N/A 


N/A 


Left 


Right 


Yes 


Left 


Right 


No 


Left 


Le ft 


Yes 


Left 


Left 


No 


Right 


Right 


Yes 


Right 


Right 


No 


Right 


Left 


Yes 


Right 


Left 


No 




CNTACT 




1 contact surface 


N/A 


N/A 


Left 


Right 


Yes 


Left 


Right 


No 


Left 


Left 


Yes 


Left 


Left 


No 


Right 


Right 


Yes 


Right 


Right 


No 


Right 


Left 


Yes 


Right 


Left 


No 



contact surface 
at t+dt 



contact surface 
at t 







Figure 3,4 Schematic of SHOCK = 100 and CONTACT = 321 
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Figure 3.5 Schematic of SHOCK = 221 and CNTACT = 222 
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the code examines the shock and contact surface condition, 
along with relative location of the shock to a contact 
surface if both exist near the node, if the flow is super- 
sonic or subsonic at the node, and the direction the flow is 
traveling. The appropriate algorithm is determined and 
called in the form of a condition subroutine. 

The "SWEEP" routine is designed to handle all 
possible shock and contact surface combinations. Because 
the code does not solve the situation where the shock and 
contact surface intersect, as illustrated in Fig. 3.6, then 
any combinations after the shock and contact surface have 
crossed call condition routines that currently contain only 
messages. The code was intentionally prepared to be easily 
extended to treat the intersection of the shock and the 
contact surface, and beyond. 



shock and 

contact surface contact surface shock 

at t .. ► intersect at t+dt ^ at t 




Figure 3.6 Schematic of SHOCK = 231 and CNTACT = 322 
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5 . The "Condition" Subroutines 

The ten subroutines that are called by "SWEEP” to 
calculate the interpolated values of the extended Riemann 
variables, R and Q, plus S, 3q/8s, and 3A/9 s are CONDI, 
C0ND2, C0ND3, C0ND4 , C0ND5 , C0ND6 , C0ND7 , C0ND7S, C0ND7N, 
and C0ND8. The procedure used in CONDI, C0ND2 , C0ND3 , and 
C0ND5 is that of Salacka [Ref. 2:pp. 36-37] modified to 
account for contact surfaces. Appendix C, Fig. C.5 is a 
general flowchart for these four routines. Essentially, an 
initial estimate is made of the characteristic slopes, X, by 
enforcing the principle of domain of dependence. The 
assumption is made that the slopes are linear, and that the 
characteristics pass through point B, as defined in Fig. 
3.7. Then q and A are computed at point A, which in turn 
yields a second estimate of the characteristic slope, A. 

The two slopes for each characteristic are compared. 
If they are not equal to within an acceptable error, the new 
estimate of A is used to repeat the process until 
convergence is achieved. All three characteristics are 
handled simultaneously. The value of 

As = A'3t (3.2) 

is then determined for each characteristic q+A, q, and q-A. 
This allows linear interpolation of Q, R, and S at point A. 
The spatial derivatives, 3Q/3s and 3V 3s are computed from 
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t+dt ◄ ► B 




Figure 3.7 Computational Method for CONDI Routine 

the values of q and A used in the estimate of the initial 
characteristic slopes. 

The "CONDI” subroutine is used when there are no 
contact surfaces or shocks within the interval on either 
side of the node, and flow is subsonic. The q+A character- 
istic uses forward differencing, while the q-A characteristic 
uses a backward differencing scheme to keep the initial 
domain of dependence of the numerical scheme outside the 
physical domain of dependence. 

The "C0ND2" subroutine is a backward differencing 
algorithm used in situations when one or more discontinui- 
ties are in the right interval, or if flow is supersonic to 
the right. Fig. 3.8 illustrates the computational method 
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t+dt B discontinuity 




Figure 3.8 Computational Method for C0ND2 Routine 

for ''C0ND2”. The characteristic slopes of q+A and q are 
handled as in "CONDl'l# but the q- A characteristic slope is 
determined by a baclward rather than a forward scheme. To 
interpolate between node i and i+1 would be incorrect 
because of the discontinuity in the values between them. 
The value of parameters in. the right inter'/al are interpo- 
lated from values in the left interval by assuming the 
derivatives do not change between the intervals up to the 
discontinuities (i.e., a shock and/or a contact surface). 

The ''C0ND3" -subroutine is a forward differencing 
algorithm used in situations when one or more discontinui- 
ties are in the left interval, or if flow is supersonic 
traveling to the left. Figure 3.9 shows the computational 
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ds 



Figure 3.9 Computational Method for C0ND3 Routine 

method for ''C0ND3”.- The method is similar to that in C0ND2 
but the q4A characteristic slope has to be interpolated from 
values of node i instead of node i-1 using a forward 
difference scheme. The characteristic slope for q-A and q 
are estimated by forward difference schemes using values of 
node i+1 and i respectively. 

For conditions, such as that illustrated in Fig. 
3.10, where a discontinuity is in the interval opposite the 
direction of supersonic flow then C0ND5 subroutine is 
called. A mix of C0ND2 and C0ND3 is used. For flow to the 
right the method of C0ND2 is implemented, while flow to the 
left is the same as for C0ND3 . 
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Figure 3.10 Example Condition for "C0ND5" Routine 

Subroutine "C0ND4" is called whenever a node is to 
be crossed by a discontinuity from either direction. The 
flowchart for ”C0ND4” is shown in Appendix C, Fig. C.6. Two 
integer vector arrays, LNODE and RHODE, are used to store 
the following information on the node being crossed. 

A) the node number, I 

B) the value of SHOCK 

C) the value of CNTACT 

D) the current value of the integer time step, J 

The information is used in "CORRCT" to update the node to 
the next time level. In the current code the maximum number 
of nodes that can be crossed during any time step is two. 
Figure 3.11 shows an example of this situation, and defines 
the assignment of values. Because the nodes are swept from 
left to right, the left-most node is assigned to LNODE, 
while RNODE applies to any node crossed to the right of the 
LNODE node. The values of A~w are set to zero, so 
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contact surface shock 




Figure 3.11 Example of "C0ND4" Routine Situation 

that these nodes are skipped when updating occurs in 
"SWEEP”. The main program interprets the value of J in 
LNODE and RNODE to determine if zero, one, or two nodes need 
to be corrected in "CORRCT". If zero nodes are crossed 
during a time step then "CORRCT" is not called. 

The "C0ND6" and "C0ND8" subroutines are designed for 
future extensions of the entire code. Until the Riemann 
problem illustrated in Fig. 3.6 is coded, then any situation 
after the shock and contact surface intersect cannot be 
computed. However, "SWEEP" is already coded to call "C0ND6" 
for the situation illustrated in Fig. 3.12, or its mirror 
image. "C0ND8" would be called for all the other situations 
after the shock and contact surface interaction has taken 
place. Currently, both subroutines output messages on the 
screen describing the situation, and set HALT equal to 1 to 
terminate the program. _ 
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Figure 3.12 Example of "C0ND6" Situation 

The ''CONDV, ''C0ND7N", and "C0ND7S" subroutines are 
all related to each other. They are intended to contain 
code that would facilitate the shock and contact surface 
intersection and solve the resulting Riemann problem. The 
"C0ND7" routine is called when the contact surface and shock 
moving in opposite directions would cross in a time step 
within an interval. The routine would calculate when and 
where within the time and space domain the intersection 
would occur based on the known velocity of each discontinu- 
ity and their current locations. This new time could then 
be used to rerun "SWEEP" with the shock and contact surface 
exactly on top of each other at the end of the new time 
step. The "C0ND7N" routine is for the same situation, but 
when a node is crossed by either discontinuity during the 
same time step. The same procedure is done, with the 
requirement to ensure the node crossed is properly updated 
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by •'CORRCT” during the new time step. The "COND7S" routine 
is called when a shock and contact surface are located at 
the same point at a time other than zero. The subroutine 
would contain the code to solve the Riemann problem set up 
by ''C0ND7N" or ''C0ND7'' routines. 

6. The "CORRCT" Routine 

Appendix C, Fig. C.7 shows the flowchart for 
"CORRCT". The "CORRCT" routine corrects nodes that have 
been crossed by a discontinuity (i.e., a shock or a contact 
surface) or are straddled by both. The routine takes the 
information from LNODE and RHODE arrays to determine which 
node or nodes are to be updated, and if there is any shock 
and contact surface interaction at the nodes. Then, using 
subroutines "DELTAX", "SXJUMP", "CSJUMP", "EXTRAP", 
"INTERP", and "BBDRY" the concepts from Section II are 
imposed to calculate the jump in the parameters R, Q, S, A 
and q at the node in question to the new time level. 

7. The "BONDRY" Routine 

The "BONDRY" routine computes the new values of the 
parameters Q, R, and S at the boundary nodes for time step 
t. The concepts of Section II for an open or closed 

boundary are coded as shown in the flowchart of "BONDRY" 
given in Appendix C, Fig. C.8. The routine uses the 
"CONDI", "C0ND2", or "C0ND3" routines to calculate the 
correct w, and then solves for w the same as in "SWEEP". 
If a shock crosses an open boundary node during a time step. 
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then "SKJUM?" routine is used to calculate the values behind 



the shock at the node. 

8. The "SRFLCT” Routine 

The "SRFLCT" subroutine is called to calculate the 
values of Q, R, S, q and A at the boundary node when a shock 
reflects from a solid boundary. Appendix C, Fig. C.9 shows 
the flowchart of "SRFLCT" which codes the analysis in 
Section II. D. 2 for closed boundary shock reflection. The 
time for the shock to reach the solid boundary, 3t^/ is 
computed. Then the excess time in the time step, 3t, is 
given by 



3tex = 9t - 3t^ (3.3) 

After the new shock velocity, Vg, is calculated, then the new 
location of the reflected shock is known from multiplying Vg 
by 9 tgj^ and adding or subtracting from the boundary 
location. 

9. The "DELTAX" Routine 

The "DELTAX" subroutine is used in "CORRCT" and 
"BONDRY" to calculate the location within an interval of a 
discontinuity in terms of x. Figure 3.13 defines the 
various displacements which are calculated. 

10. The "SKJUMP" Routine 

The "SKJUMP" subroutine computes the conditions 
downstream of a stationary normal shock as described in 
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Figure 3.13 Discontinuity Location within an Interval 

Section II. D. 2. The "CORRCT", "TRAK", and "BONDRY" 
subroutines call "SKTUMP" when required, to obtain the speed 
of sound (A) , entropy (S) , and velocity (q) behind the 
shock. 

11. The ”CSJUMP» Routine 

The "CSJUMP" subroutine computes the velocity and 
speed of sound change across a contract surface as described 
in Section II. D. 3. The "CORRCT" and "TRAK" routines call 
"CSJUMP". 

12 . Numerical Support Routines 

The "EXTRAP", "INTERP", and "BBDRY" routines are 
used by "CORRCT", "BONDRY", "TRAK", and "DBURST" to extrapo- 
late or interpolate data to the face of a discontinuity or 
point. In particular, "BBDRY" is an interpolation routine 
used within an interval near a boundary. 
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13 . 'The Output Routines 



The four output subroutines used are those developed 
by Salacka [Ref. 2;pp. 39-40]. The "BORDER" and "PLOT" 

routines output plots of pressure, velocity, density, and 
modified entropy distributions versus a non-dimension- 

alized tube length at preset time intervals. "EXACT" 
creates a plot of the exact density distribution at selected 
points and compares it with the computed density 

distribution for a selected time step. The "LIST" routine 
outputs to files 8, 9, and 10 tabular listing of computed 

data. The files are sent to the user's permanent storage 
disk. 

A value of GRAPHS equal to one causes "BORDER" to 
be called in the main program. This defines plot axis, 

labels, and headings. "PLOT" is then called once at time 
zero, and at every time step set by SKIP. 

If GRAPHS is set equal to two, then "EXACT" is 

called every SKIP time steps to plot six exact density 

values and the computed density distribution. The six exact 
points are: 

A) The two boundary points 

B) The head and tail of the rarefaction wave 

C) A point just behind the contact surface 

D) A point just behind the shock. 
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The location of the exact values is based on elapsed time 
and known exact values for wave velocities entered with the 
initial conditions for the problem. 

"LIST" routine is called every SKIP time steps when 
GRAPHS is equal to zero. File 9 contains a tabular listing 
of the Riemann variables, modified entropy, pressure, 
density and velocity distributions, elapsed time, discon- 
tinuity velocity and location. File 10 is a tabular listing 
of the location of the shock, contact surface, and the head 
and tail of the expansion wave computed at every SKIP time 
steps. File 8 contains the exact values of the location of 
the shock, contact surface, and the head and tail of the 
expansion wave. 
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IV. RESULTS 



Four test cases were run using the E1DV2 code. The 
purpose was to verify the ability of the code to determine 
the unsteady flow process and correctly simulate varying 
boundary conditions and flow directions in the Riemann shock 
tube problem. The shock tube is illustrated in Fig. 4.1 
with the high pressure on the left. The tube is divided 
into two sections by a diaphragm. When the diaphragm is 
burst, the pressure equalizes through a shock wave traveling 
into the expansion chamber, and an expansion or rarefaction 
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P 

low 



time = 0 




Figure 4 . 1 Shock Tube at t = 0. and t = t 
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wave traveling into the compression chamber. A contact 
discontinuity is behind the shock wave, traveling at the 
particle velocity [Ref. 6:p. 30]. 

A. TEST CASE 1 

Test Case 1 was designed to demonstrate shock reflection 
and expansion wave reflection at solid boundaries. Test 
Case 1 had the following initial and boundary conditions: 

1) Pressure and density ratios equal to five, with high 
pressure on the left 

2) Temperature ratio equal to unity 

3) Both boundaries were closed 

4) Diaphragm was located at x = 0.5 

5) Computational mesh had 101 nodes 

6) Maximum time step, JSTOP, = 109, with SKIP = 18. 

Plots of the results obtained for the pressure, density, 
velocity, and modified entropy distributions are shown in 
Fig. 4.2 and Fig. 4.3. Fig. 4.2 shows the computation up to 
time step 55, while Fig. 4.3 takes the computation from time 
step 56 to termination. The output of "EXACT", a plot of 
the exact density compared with the computed density 
distribution is shown in Fig. 4.4. Fig. 4.5 illustrates the 
spatial location versus time for exact and computed shock, 
contact surface, and head and tail of the expansion wave. 
Data for Fig. 4.5 were taken from file 10 and file 3. Fig. 
4.6 shows plots of pressure, density, modified entropy, and 
velocity distributions for the same initial and boundary 



77 



PRESSURE VELOCITY 

0.0 2.0 4.0 6.0 - 2.0 - 1.0 0.0 



SHOCK TUBE RESULTS 

FIRST ORDER N -101 

DENSITY RATIO - 5.0 TEMP RATIO - 1 
PRESSURE RATIO - 5.0 



o in 




i i I I 1 i i i i 

0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 




Figure 4.2 Test Case 1 (J = 1 to J = 55) 
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Figure 4.3 Test Case 1 (J = 56 to J = 109) 
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Figure 4.4 Exact and Computed Density Distributions 
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Figure 4.5 Location of Discontinuities versus Time (Test Case 1) 
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Figure 4.6 Test Case 1, High Pressure on Right 
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conditions but with high pressure on the right. Figure 4.6 
includes data corresponding to those in both Fig. 4.2 and 
Fig. 4.3 for the high pressure initially on the left. 

B. TEST CASE 2 

Test Case 2 was to demonstrate an open boundary, 
constant pressure, expansion wave interaction. The test 
case was run with the following initial and boundary 
conditions ; 

1) Pressure and density ratios equal to 5.0, with high 
pressure on the left 

2) Temperature ratio equal to unity 

3) The left boundary was open, with a constant pressure 
ratio across the boundary of 4/5 

4) The right boundary was open, with an adjustable pres- 
sure ratio that in effecr extended the length of the 
tube 

5) The diaphragm was located at x = 0.74 

6) Computational mesh had 51 nodes 

7) Maximum time step, JSTOP, = 37, SKIP = 9. 

Plots of the results for the pressure, density, velocity, 
and modified entropy distributions are shown in Fig. 4.7. 
The test case was rerun with high pressure on the right. 
Figure 4.8 shows the results for pressure, density, and 
temperature ratios the same. The diaphragm was then at x = 
0.26. Boundary conditions were reversed. The same computa- 
tional mesh, and output control were used. 
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Figure 4.7 Test Case 2 , High Pressure on Left 
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Figure 4.8 Test Case 2, High Pressure on Right 
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C. TEST CASE 3 

Test Case 3 was designed to demonstrate a shock exiting 
an open boundary with constant pressure maintained at the 
boundary itself. The following initial and boundary 
conditions were set; 

1) Pressure and density ratios equal to 5.5, with high 
pressure on the left 

2) Temperature ratio equal to unity 

3) Left boundary was closed 

4) Right boundary was open with a constant pressure ratio 
across the boundary of unity 

5) Diaphragm was located at x = 0.5 

6) Computational mesh had 101 nodes 

7) Maximum time step, JSTOP, = 73, with SKIP = 18. 

Figure 4.9 shows plots of the results obtained for pressure, 
density, modified entropy, and velocity distributions. 
Similarly, Fig. 4.10 shows results obtained by putting the 
high pressure on the right, reversing the boundary condi- 
tions, and holding everything else the same. 

D. TEST CASE 4 

Test Case 4 was designed to demonstrate a lower pressure 
ratio, with shock wave reflection. The initial and boundary 
conditions were set as follows: 

1) Pressure and density ratios equal to 3.2, with high 
pressure on the left 

2) Temperature ratio equal to unity 

3) Both boundaries were closed 
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Figure 4.9 Test Case 3, High Pressure on Left 
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Figure 4.10 Test Case 3, High . Pressure on Right 
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4) Diaphragm was located at x = 0.5 

5) Computational mesh had 901 nodes 

6) Maximum time step, JSTOP, = 601, with SKIP = 150. 

Plots of the results obtained for the pressure, density, 
modified entropy, and velocity distributions are given in 
Fig. 4.11. Figure 4.12 shows the results obtained with the 
high pressure to the right, SKIP = 200, and holding all 

other parameters constant. 
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Figure 4.11 Test Case 4, High Pressure on Left 
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Figure 4.12 Test Case 4, High Pressure on Right 
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V. DISCUSSION 



A. RESULTS OF TEST CASES 

Overall, in the four specific test cases which were run, 
the expected qualitative flow behavior was produced by the 
code. Quantitatively, either experimental data or further 
exact solutions are needed to fully validate the computa- 
tions. The results of each test case will be discussed 
separately. 

Test Case 1 results in Fig. 4.2 show the formation of a 
well-defined shock wave traveling to the right, with the 
contact surface crisply defined following along behind. The 
entropy drops sharply across the shock and remains constant 
to the contact surface and then jumps to a steady, constant 
value to the left boundary as expected. Pressure and 
velocity remain constant through the contact surface. 
Velocity is positive since flow is to the right. The 
expansion wave is smeared from head to tail over the correct 
range. Figure 4.3 is a continuation of the flow problem 
with a sequential display of the results. The entropy 
continues to drop as the reflected shock travels back toward 
the contact surface. Notice that the velocity behind the 
shock is zero. At the left boundary as the head of the 
expansion wave reflects, pressure and density continues to 
drop while velocity is zero at the boundary. The density 
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distribution shows that the contact surface and shock are 
about to cross. The program terminated and issued a message 
before the next output call was made, when the code 
determined that the contact surface and shock would cross. 
The exact density distribution compared with the computed 
density distribution in Fig. 4.4 clearly shows the shock and 
contact surface are modeled very accurately with a medium 
mesh. The expansion process, consisting of infinitesimal 
changes propagating along the characteristics is fairly well 
modeled. Figure 4.5 shows that tracking of the contact 
surface and shock wave, in comparison with exact locations at 
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slightly different from those of an exact Riemann solver 
code [Ref. 8] . The exact code is based on a method given in 
[Ref. 9:pp. 181-191]. However, an examination of the 

tabulated output showed that the solution method predicted 
changes in the gradients to occur at similar locations to 
those computed using the q+A and q-A estimates. Consequent- 
ly, while there may be some inaccuracy in the numerical 
solution, the method of tracking the head and tail of the 
expansion wave is acceptable. 

The ability to reproduce the same conditions but with 
flow to the left is demonstrated clearly in Fig. 4.6. The 
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entropy, density, velocity, and pressure distributions are 
mirror-images of those for the case of high pressure on the 
left. 

Test Case 2 results in Fig. 4.7 show that at time zero 
with a pressure ratio at the boundary of 4 , which is lower 
than the preset value inside the tube of 5, an expansion 
wave traveling inwards is generated. Outflow is seen in the 
negative velocity distribution developing near the left 
boundary, where negative velocity implies flow traveling to 
the left. The contact surface and shock are clearly formed 
to the right of the off-center diaphragm, and travel to the 
right. The expansion wave generated at the diaphragm 
travels to the left and intersects with the expansion wave 
generated at the left boundary. The pressure distribution 
shows the propagation of these waves as the pressure drops 
behind the head of each wave. With conditions reversed, so 
that the high pressure is to the right of the diaphragm, the 
results in Fig. 4.8 are a mirror- image with velocity 
negative for flow to the left, and positive for the outflow 
to the right. 

Test Case 3 results in Fig. 4.9 demonstrate that the 
shock wave meeting an open boundary where the pressure is 
held constant, is correctly computed to exit the tube. The 
density distribution shows the jump increase across the 
shock, then another jump increase across the contact 
surface. Then density plunges down as an expansion wave 
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travels back into the tube when the shock exits. The 
pressure ratio at the open right boundary was held at unity, 
as seen in the pressure distribution which, at the boundary 
behaves similarly to the density distribution. The pressure 
at the boundary would be required to increase as the shock 
passed by, however the enforcement of the constant pressure 
locally causes an expansion wave to form traveling to the 
left. The entropy remains constant, at the value behind the 
shock, from the boundary back to the contact surface. There 
is a jump in entropy across the contact surface. The 
velocity distribution shows that as the shock travels to the 
right, the velocity jumps up. When the expansion wave 
generated at the right boundary travels inward the velocity 
continues to increase in magnitude while flowing to the 
right. The expansion wave generated at the diaphragm can be 
seen traveling to the left in the plots of pressure, 
density, and velocity. Reversing the situation and 
computing conditions for flow to the left, in Fig. 4.10, 
results in a mirror- image in the pressure, density, and 
entropy distributions. The velocity becomes negative since 
flow is to the left. 

Test Case 4 is a repeat of the first test case but with 
a pressure ratio of 3.2 and a fine computational mesh. In 
Fig. 4.11 the density steps up across the shock and contact 
surface as they travel to the right. When the shock 
reflects the density jumps again. The velocity distribution 
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shews that the velocity drops to zero behind the reflected 
shock, after it had jumped up across the shock when it was 
traveling to the right. The instability seen in the 
velocity and pressure distributions near the midpoint 
occurred during the first 100 time steps. A numerical 
instability appears to generate at the diaphragm at time 
zero for pressure ratios less than about 3.0 which can be 
severe. For the conditions in this test case the transient 
instability damped out after the first 100 time steps. The 
shock, contact surface, and expansion wave are nevertheless 
seen to be accurately modeled. Figure 4.12 demonstrates 
that reversed conditions result in mirror-images of the 
computed conditions. 

B. CURRENT LIMITATIONS OF THE CODE E1DV2 

The current limitations of the E1DV2 code for modeling 
quasi-one-dimensional inviscid flow are in two categories. 
First, there are numerical limitations in obtaining 
solutions with the present code for different initial and 
boundary conditions. Second, the present coding limits what 
flow situations can be treated. 

The numerical instability which occurs at low pressure 
ratios during the first few time steps can be severe. The 
assumption in the current code is that the shock forms 
immediately, which at high pressure ratios does not pose any 
problems. The mathematical modeling of the bursting 
diaphragm in subroutine "DBURST" adequately reduces the 
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transient instability for high pressure ratios, but not for 
low pressure ratios. 

A second rather different numerical limitation was 
identified in numerically detecting the location of the head 
and tail of the expansion process. Because the expansion 
process is not a sharply defined front, fixing the precise 
locations of the head and tail waves numerically at a given 
time is difficult. However, the method currently used does 
accurately track the computed propagation of the wave along 
the characteristics. The characteristics are modeled 
currently as straight lines. A higher order curve would 
describe them more accurately. 

The current code is limited to tracking two discontinui- 
ties and the expansion wave (i.e., a shock and a contact 
surface) . Thus any situation which would generate another 
discontinuity can not be computed. The code can determine 
when this will occur, and will issue a message explaining 
the situation before terminating. Thus the shock colliding 
with a contact surface will currently terminate the program. 
A boundary pressure that is higher than the pressure inside 
the tube would also create a compression wave or possibly a 
shock wave. This condition will also terminate the program. 
Simulation of the process of a shock forming over an 
extended distance and time as compression waves pile on top 
of each other and- strengthen, would also require additional 
code. 
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Finally, the variation of gamma that occurs across the 
contact surface in a wave rotor cannot be handled currently 
since E1DV2 assumes a single constant gamma throughout. 
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VI. CONCLUSIONS 



Towards the development of a one-dimensional code for 
wave rotor applications, methods for tracking and correcting 
conditions across a contact discontinuity, applying open and 
closed-end boundary conditions, accounting for shock wave 
and contact surface interaction were devised and were 
presented here in detail. The EULERl Fortran Code [Ref. 2] 
was revised to become the E1DV2 code with the following 
additional capabilities: 

1) tracking of the contact surface and expansion wave 

2) imposing high pressure initially on the right side of 
the diaphragm 

3) correct jump conditions across the contact surface 

4) allow open boundary conditions with constant pressure 
specified and allow exiting of shock and expansion 
waves 

5) allow closed boundary conditions and model shock and 
expansion wave reflections 

6) improved numerical accuracy from time zero to the 
maximum time step. 

The code was tested on the shock tube problem under four 
different initial and boundary conditions with excellent 
results . 

The following extensions are recommended to make the 
code suitable for wave rotor applications: 

1) Solve the Riemann problem at the moment when the shock 
and contact surface intersect 
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2) Add additional code to track more than two discontinu- 
ities and one expansion wave 

3) Incorporate a variable value of gamma into the code 

4) Update the characteristic curve approximation from 
linear to a higher order polynomial 

5) Add code necessary to handle open boundary conditions 
with inflow 

6) Improve on the numerical computation at time zero for 
low pressure ratios across the diaphragm 

7) Enable boundary conditions to be variable during 
program execution to allow wave rotor cycle design 

8) Compare computations using the code with experimental 
data where available 

9) Extend the code to two dimensions. 

These extensions may require various degrees of effort. 
However, the ability of the QAZID solution method to be 
extended rather simply to describe viscous multi-dimensional 
flow suggests that such efforts would be justified. 
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APPENDIX A 



DERIVATIONS OF EQUATIONS 



A. LIST OF VARIABLES 



A 

Cv 

e 

h 

P 

Q 

qR 

q 

R 

Rg 

s 

T 

t 

u 

V 
W 

w 

p 

Y 



Speed of sound 

Specific heat at constant pressure 
Specific heat at constant volume 
specific internal energy 
Specific enthalpy 
Static pressure 
Modified Riemann variable 
Reversible heat transferred 
Velocity magnitude 
Modified Riemann variable 
Gas constant 

Modified entropy (non-dimensional where S = S-Rq) 

Static temperature 

Time 

Velocity relative to a standing shock wave 
Specific volume 

Incoming Mach Number relative to a stationary 
shock wave 

Work 

Density 

Ratio of specific heats 
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B. DERIVATION OF SHOCK JUMP EQUATION FOR HIGH PRESSURE ON 

THE RIGHT 

The analytical equation relating the non-dimensionalized 
extended Riemann variable, R, change through a normal shock 
is derived below similar to that for the Q variable change 
when high pressure is on the left [Ref. Z :Appendix A]. 

Figure A.l shows a shock moving to the left with 
velocity Vg. Subscript A is always associated with 
parameters to the left of any discontinuity, and B with 
those on the right. To enable normal shock relations to be 
used the system requires a velocity in the opposite direc- 
tion but of equal magnitude be imposed upon it. The 
coordinate system was defined such that vectors, such as 
velocity, directed to the right are positive while those to 
the left are negative. 
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Figure A. 1 Shock Wave with High Pressure on Right 
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since relative incoming Mach Number, w, relates the 
velocity downstream, or that region into which the shock is 
traveling, to the speed of sound in that region [Ref- 10:pp. 
114-154] then 



w 



S 







(Al) 



w is a positive value because though negative, is small- 
er in magnitude than the equally negative shock velocity, 
Vg. The speed of sound. A, is positive by definition. 

The appropriate extended Riemann variable in this case 
is R, defined as 



R = q - As (A2) 

Since q is negative, this can be rewritten as 

R = -(|q| + AS) (A3) 

to emphasize that R is the Riemann variable associated with 
the lesser change across the shock [Ref. 3:pp. 18-21]. We 

adopt the pattern of non-dimensionalizing velocity by the 
speed of sound, pressure and density by corresponding values 
downstream of the shock. Using the entropy S, defined by 
Verhoff [Ref. l:p. 2] as 
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(A4) 



then 
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(A5) 



The ratios of pressure, density, and sonic velocity across a 
normal shock from Zucker [Ref. 7:p. 151] and Shapiro [Ref. 
10 :p. 113] in terms of Mach Number are 



*E 






, 2 y ,.,2 _ ri 

Y+1 Y+1 



^ = (Y+l)w^ 

Pa (y-1M^ +2 



(Y+1)w 



[2 (y- 1) [1 -l]-3 



1/2 



(A6) 



(A7) 



(A8) 



The first term on the right side of equation (A5) can be ex- 
pressed in terms of w by applying simple continuity in mass. 
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P^UAArea^ = PBUBAreag 
Taking the areas as egua], this becomes 



(A9) 
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Combining equations 
equation (A'5) gives 



(A6) , (A7) , 



(A8) , and (A13) into 
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C. DERIVATION OF CONTACT SURFACE JUMP EQUATION 

The analytical expression for the ratio of sonic 
velocity across a constant surface is derived / as follows; 
The first law of thermodynamics is stated as 



■ 




incremental 




work done ' 


change of 




amount of 




on the 


internal energy 


= 


heat added 


+ 


system by 


■ _ 




to system 




surrounding 



or, in differential form 



de = 3qR + 3 W (A15) 

Internal energy is defined as 



e = h - pv (A16) 



then 
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de 



dh - pdv - vdp 



(A17) 



The incremental work is given by 

3W = pdv (A18) 

The heat addition is determined from the definition of modi- 
fied entropy. 



dS = “ — 



1 '^R 



Y T 



dqR = - yTdS 



(A19) 



Substituting equations (A17) , (A18) , and (A19) into equation 

(A15) yields 

dh - pdv - vdp = -yTdS - pdv 
Canceling like terms, and rearranging gives 



dp = 1 TdS + (-)dh 

V V 



(A20) 



Enthalpy, h, is defined as 



h = CpT 



(A21) 
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For a perfect gas the sonic velocity is given by 



A = 



thus 



dA = RGT)^/^(dT/T) (A22) 

Substituting equation (A21) and then (A22) into equation 
(A20) , equation (A.20) becomes 

dp = ^TdT + (i)d(CpT) 

= iTdS + (i) CpdT 

= iTdS + (i)CpT(^) (A23) 

For an ideal gas, 

Y = Cp/Cy 



and 



Cp — Cy + Rq 

Combining these two equations, and rearranging 
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(A24) 




Substituting equation (A24) into equation (A23) , and 
observing that pressure remains constant through the contact 
surface 



-yds 

V ' 



— (R f ) ) 
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Eliminating T/v, equation (A25) becomes 



(A25) 
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Then S/Rq becomes a non-dimensionalized entropy. Using the 
notation of S = S/Rq now for the non-dimensionalized form, 
so that by integrating both sides 
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which can be written as 
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exp 






(A27) 



Definitions used in this development, excluding modified 
entropy were from [Ref. ll:pp. 87-125]. 
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D. DERIVATION OF EXTENDED RIEMANN VARIABLE CHANGE ACROSS A 
CONTACT SURFACE 

An analytical expression for the change in the non- 
dimensionalized extended Riemann variable, Q, across a 
contact surface traveling right is derived here. Figure A. 2 
depicts a contact surface moving right with velocity q. 



Entropy 




Figure A. 2 Contact Surface Traveling Right 
Subscript Notation 

Values of parameters to the left of the interface are 
denoted by subscript A', and those to the right with sub- 
script B'. As illustrated, parameters to the left of the 
shock have subscript A, and to the right subscript B. All 
velocities are non-dimensionalized by sonic velocity 
immediately downstream of the discontinuity. Thus using the 
definition of the extended Riemann variable, 
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Q 



q + As 



(A28) 



then 







“V 






(A29) 



Since velocity is constant through a contact surface, so 
that 



^A • “ <3b ' 



then equation (A29) becomes 
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Using equation (A27) , this can be written 
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= exp 
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(A30) 



Multiply each side by A_,/A„, and since A_ , = A then 

o B BA 



^A'”^B' 






(■^^(Sg.-s .) 
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(A31) 
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For a contact surface traveling to the left, as 
illustrated in Fig. A. 3, the derivation is entirely similar but: 




Figure A. 3 Contact Surface Traveling Left 
Subscript Notation 



using the R, extended Riemann variable with the result 






= [exp 
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E. DERIVATION OF OPEN BOUNDARY CONDITIONS WITH CONSTANT 
PRESSURE 

The equation for density and temperature ratios in terms 
of pressure and entropy are derived below. The definition 
of modified entropy in non-dimensional form is 



S 



2 
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p" 



(A33) 
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where P and p are non-dinensionalized ratios. Then 



(S - (-Y(Y-D) = lJ^(P/P^) 



or 



exp ^ ^ = P/P^ 



Rearranging, 






t(S --4-)(-Y(Y-1))3 
[P [exp ' ] ] 



thus 



1/y (S-— j-)(-Y(Y-l)) -Vy 

p = P'^'^(exp ^ ) 



, (s-X)(-y(y-D) -Vy 

p = {|(exp Y-i )} 



using the perfect gas law» • 
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) 



multiply each side by T and divide by p ^ then 



(A34) 
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T 



thus 



Y (S-;^)(-Y(Y-l)) 

= p VP (exp ^ ) 



Y-1 (S-;^)(Y(1-Y)) 

‘ = "^[exp ^ ] 



(A35) 
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APPENDIX B 



OPERATION OF E1DV2 ON THE NPS VM/CMS SYSTEM 
A. TERMINAL REQUIREMENTS 

Terminal requirements depend on the desired output. Any 
terminal, such as the IBM 3278, connected into the VM/CMS 
system is adequate for a tabular listing or Disspla Metafile 
of data. The Disspla Metafile stores graphical data for 
display using DISSPOP commands. If output of the plots on a 
monitor screen is desired, an IBM 3277-TEK618 dual screen 
terminal must be used. Table B.I lists terminal require- 
ments for different outputs. 



TABLE B.I 

TERMINAL AND OUTPUT 



OUTPUT 



TERMINAL 



Tabular listing of data Any terminal 

Graphical data in format 

for VERSATEC printing Any terminal 



Graphical data plotted 

on screen IBM 3277-TEK618 



B. PROCEDURE 

1. Editing Variables and Program Setup 

To change the initial and boundary conditions, 
graphical output, and grid size to run the program under 
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different conditions^ the following must be done. First, 
pick the appropriate terminal from Table B.I and log on. 
Enter the editor by issuing the command 

X E1DV2 FORTRAN A 

Table B.II lists the code variables that will require 
editing and their meaning. In addition, to these variables 
if graphical plots are to be outputed then enter the 
"BORDER” and "EXACT" subroutines and edit the lines: 

"FIRST ORDER N = ???" 

"DENSITY RATIO = ??? TEMP RATIO = ???" 

"PRESSURE RATIO = ???" 

To issue the output graphs to the screen comment out the 
lines "COMPRS", and use 

CALL TEK618 

Otherwise, if a Disspla Metafile of the graphics plot is 
desired, comment out the "TEK618" line, and use 

CALL COMPRS 

These lines are in the "Set up graphics plot of variable" 
loop in the main program. 

To store these changes, hit the enter key, and type 
"file". Select the enter key once more. The program is now 
ready to be compiled and executed. 

2 . Commands 

To compile and execute the E1DV2 code on the VM/CMS 
system perform the following commands exactly as typed: 
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TABLE B.II 



Variable 



EDITABLE VARIABLES 
Edit 



DIMENSION 


statement 


set the dimension of the arrays 


for arrays 


A(N) 


equal to the number of nodes 


XARRAY(N) 




in the grid 


N 




set to the number of nodes in the 
grid 


GRAPHS 




set to: 0 for tabular listing of 

data 

1 for plot of density, 
entropy, pressure and 
velocity distributions 

2 for plot of exact density 
distribution compared to 
computed density 
distribution 


SKIP 




set to number of time steps between 
calls to output routines 


JSTOP 




set to maximum number of time steps 


TRI 




set to initial temperature ratio 
across diaphragm 


PRI 




set to initial pressure ratio across 
diaphragm 


DRI 




set to initial density ratio across 
diaphragm 


QLI 




set to initial velocity left of 
diaphragm 


QRI 




set to initial velocity right of 
diaphragm 


LBDPRI 




set to initial pressure ratio 
across left boundary 


LBDTRI 




set to initial temperature ratio 
across left boundary 


LBDORI 




set to initial density ratio across 



left boundary 
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TABLE B.II (CONTINUED) 



RBDDRI 

RBDPRI 

RBDTRI 

G 



set to initial density ratio across 
right boundary 

set to initial pressure ratio across 
right boundary 

set to initial temperature ratio 
across right boundary 

set to desired value of gamma 



EE 



LWPRES 



LBNDRY 

RBNDRY 

LBDPRS 

RBDPRS 

XINIT 

VHEAD 

VTAIL 

VCDE 

VSE 



set to desired error tolerance for 
calculation of characteristic slope 
(i.e. , O.lD-8) 



set to: 2 for low pressure on right 

side of diaphragm 

3 for low pressure on left 
side 

set to: 1 for closed boundary 

0 for open boundary 

set to: 0 for constant pressure at 

left boundary 

1 for pressure that adjusts 
at the left boundary 

set to location of diaphragm 



set to exact velocity for head of 
expansion wave 

set to exact velocity for tail of 
expansion wave 

set to exact velocity for contact 
surface 



set to exact velocity for shock 
wave 



DLCD 



set to exact density behind contact 
surface 



DLSH 



set to exact density behind shock 
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TABLE B.II (CONTINUED) 



SIGMA (1,2) 
SIGMA(2,2) 
SIGMA(3,2) 
SIGMA(4,2) 

Y 



set to diaphragm location (i.e., 
0.5D00) 



There are four cases where Y 
appears in the program 
edit as follows: 

Comment out Y = (Integer#) , use: 
first, Y = (N+l)/2 if LWPRES = 2 
second, Y = (N+3)/2 
and 

third, Y = (N-l)/2 if LWPRES = 3 
fourth, Y = (N+l)/2 
for diaphragm at 0.5 
Comment out those above, if dia- 
phragm at another node. Set 
first, Y = (Integer #) of node for 
LWPRES = 2 
second, Y = (# + 1) 
and 

third, Y = (Integer # - 1) of node 
for LWPRES = 3 
fourth, Y = (#) 
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1) Increase the virtual memory by entering 

DEFINE STORAGE IM 

2) To return to CMS environment enter 

I CMS 

3) To compile the program enter 

FORTVS E1DV2 

The screen will display messages as it compiles each 
routine and when finished a ready symbol appears. 

4) To execute the program, enter 

DISSPLA E1DV2 

The message 

...YOUR FORTRAN PROGRAM IS NOW BEING LOADED... 

. . .EXECUTION WILL SOON FOLLOW. . . 

should appear, followed by 

. . .EXECUTION BEGINS. . . 

If at a TEK618 terminal with GRAPHS equal to 1 or 2 
then the screen on the TEK618 will begin plotting the 
selected graph. A 

...press ENTER to continue... 

message will appear on the 3277 terminal. If a copy of the 
plot is desired, do so now before pressing the enter key. 
After pressing the enter key on the 3277 terminal, the plot 
will be erased and the program will terminate. Proper 
termination will result in 

END OF DISSPLAY 9.2 #### VECTORS IN 1 PLOT... 
appearing, followed by a ready message. 
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If GRAPHS was set to 0, then proper termination 
would be a ready message. The tabular listing of pressure, 
density, velocity, entropy, and Riemann variables will be in 
"FILE FT09F001." The exact location of the shock, contact 
surface, and expansion wave with elapsed time will be in 
"FILE FT08F001." The computed location of the shock, 
contact surface, and expansion wave with elapsed time will 
be in "FILE FTIOFOOI." 

A second method to compile and execute the program, 
plus provide the files with a name is to create the follow- 
ing EXEC file on the user's disk. 

FI 9 DISK FILE09 LISTING A (PERM 
FI 10 DISK FILEIO LISTING A (PERM 
FI 8 DISK FILE08 LISTING A (PERM 
FORTVS &1 

A possible file name and required file type would be 

RUN EXEC 

To compile the program and define the output files, enter 

RUN E1DV2 

After compiling is finished, and the ready message appears, 
enter 

DISSPLA E1DV2 

The program executes as outlined before. 
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APPENDIX C 



FLOWCHARTS 



Flowcharts for major routines in E1DV2 are shown. 
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Figure C.l Main Program Flowchart 
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no 





call extrap 


call bondry 


1 








call srflct 


i 





repeat PlocK. A 
for left boundary 




^ stop ^ 



Figure C.l (Continued) 
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Figure C.2 "DBURST" Subroutine Flowchart 
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Figure C. 2 ' (Continued) 
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( input data) 







update slgma(lj) 
to current time 






determine 12(1) 



sKdir-=left 




□ □ 



Figure C.3 "TRAK” Subroutine Flowchart • 
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Figure C.3 (Continued) 
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Figure C.3 (Continued) 
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Figure C.4 (Continued) 
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Figure C.4 (Continued) 



133 




Input Data 




Figure C.5 General "CONDI, 2, 3, and 5" Subroutine Flowchart 
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Figure C.6 ''C0ND4” Subroutine Flowchart 
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Figure C.7 (Continued) 
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Figure C.8 "BONDRY" Subroutine Flowchart 
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Figure C.8 (Continued) 
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Figure C.8 (Continued) 



141 









Figure C.9 "SRFLCT" Subroutine Flowchart 
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APPENDIX D 



E1DV2 FORTRAN LISTING 
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flHIK 



EULER-ID 

VERSION 2 
(E1DV2) 

THIS PROGRAM SOLVES THE EULER EQUATIONS 
EXPRESSED IN A QUA2I-0NE DIMENSIONAL 
STREAMLINE COORDINATE SYSTEM. 

AUTHOR - LT. O.T. JOHNSTON, FEB 1987 



BASED ON THE EULERl CODE BY 
T.F. SALACKA, DEC 1985 



FEATURES OF THIS VERSION 12) 

♦ ORDER OF SPATIAL OERIVITIVES - FIRST 

♦ NUMBER OF SPATIAL DIMENSIONS - ONE 

♦ DISCONTINUITIES TREATED i 

SHOCKS - YES 

CONTACT DISCONTINUITIES - YES 

EXPANSION HAVES - YES 

4 HIGH PRESSURE SIDE 

LEFT - YES 

RIGHT - YES 



»IHI« 



♦ 4444 444 4^^4 ♦♦♦♦♦ 4 4444444 4 ♦♦44444f44 ♦♦♦♦♦♦♦♦♦♦♦♦ 
4 4 

4 CONVEhmONS AND DEFINITIONS 4 

4 4 

44444444444444444444444444444444444444 ^ 4444444^4 



NON-DIMENSIONING CONVENTION- 



ALL VELOCITIES NON-DIM. BY THE SOUND SPEED ON 
THE LOW PRESSURE SIDE OF THE DIAPHRAGM. 

ALL PRESSURES, DENSITIES, AND TEMPERATURES ARE 
NON-DIM. 8Y THEIR INITIAL VALUES ON THE LOW 
PRESSURE SIDE OF THE DIAPHRAGM. 

SPACIAL DISTANCE IS NON-DIM. BY OVERALL LENGTH. 
ENTROPY IS NON-DIM. BY THE GAS CONSTANT, R. 

TIME IS NON-DIM. BY < LENGTH/SOUND SPEED). 
VELOCITIES AND DISTANCES ARE DEFINED POSITIVE 
TO THE RIGHT. 



SUBSCRIPT NOTATION- 



Z 

J 

K 



L - 



SPACIAL NODE Cl TO N FROM LEFT TO RIGHT) 
TIME LEVEL 10 IS THE INITIAL CONDITION) 
DENOTES WHICH CHARACTERISTIC WAVE IS BEING 
DEALT WITH: 

1 » Q4A 2 ■ Q 5 ■ Q-A 

DENOTES WHICH TYPE OF DISCONTINUITY IS 
BEING DEALT WITH: 

1 a SHOCK 2 • CONTACT DISCONTINUITY 

3 a HEAD OF EXPANSION HAVE 

4 • TAIL OF EXPANSION HAVE 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



A - SPEED OF SOUND 

»A - DENOTES THE VALUE OF A VARIABLE AT THE NODE 
TO THE LEFT OF A DISCONTINUITY. * CAN 
BE ANY VARIABLE NAME. 

AR - THE RATIO OF SOUND SPEED ACROSS A 

SHOCK, A/B( SHOCK MOVING RIGHT ) >B/A( SHOCK MOVING LEFT) 

- DENOTES THE VALUE OF A VARIABLE AT THE NODE 
TO THE RIGHT OF A DISCONTINUITY. 

BDRY - 3 DENOTES LEFT BOUND ARY, 2 THE RIGHT BOUNDARY 

COUNT - COUNTER FOR GRAPHICS ROUTINES 

DARRAY - ARRAY OF DENSITY FOR PLOTTING 

CELT - TIME STEP 

DLCD - DENSITY BEHIND THE CONTACT 

DISCONTINUITY IN THE EXACT SOLUTION. 

DLSH - DENSITY BEHIND THE SHOCK IN THE 
EXACT SOLUTION. 

DQ - THE JUMP IN VELOCITY ACROSS THE SHOCK 

DIVIDED BY THE SOUND SPEED AT B( RIGHT) OR A(LEFT) 

DRI - INITIAL DENSITY RATIO ACROSS THE SHOCK 

EE - DESIRED PRECISION FOR CHARACTERISTIC CALCULATIONS 

G - GAMMA (RATIO OF SPECIFIC HEATS) 

GRAPHS - FOR GRAPHICAL OUTPUT, 0=NONE (TABULAR) 
l=PLOT$ ALL VARIABLES 
2=COMPARES DENSITY WITH EXACT SOLUTION 
G1 - 1/(G-1) 

G2 - 2/(G-l) 

H - 1/(N-1) 

HALT - TERMINATES PROGRAM IF 1,SET BY CONDITIONS NOT CODED 
12 - NUMBER OF THE NODE TO THE RIGHT OF A 

DISCONTINUITY. 

JSTOP - NUMBER OF TIME LEVELS TO BE CALCULATED 
LBDDR - LEFT BOUNDARY DENSITY RATIO 
LBDDRI - LEFT BOUNDARY DENSITY RATIO AT TIME ZERO 
LBDPR - LEFT BOUTJDARY PRESSURE RATIO 
LBDPRI - LEFT BOUNDARY PRESSURE RATIO AT TIME ZERO 
LBDPRS - VALUE OF 0 DENOTES CONSTANT PRESSURE AT LEFT BOUNDARY 
,1 DENOTES ADJUSTABLE PRESSURE AT THE LEFT BOUNDARY 
LBDTR - LEFT BOUNDARY TEMPERATURE RATIO 
LBDTRI - LEFT BOUNDARY TEMPERATURE RATIO AT TIME ZERO 
LBNDRY - DENOTES LEFT BOUNDARY CONDITION, OPEN OR CLOSED 
LNOOE - ARRAY OF LEFT MOST NODE TO BE CORRECTED IN CORRCT 
LWPRES - DENOTES WHICH SIDE OF DIAPHRAGM HAS LOW PRESSURE 
N - NUMBER OF SPACIAL NODES (ODD NUMBER) 

ND - DOUBLE PRECISION VALUE OF N 

STORED VALUES OF FOR THE NEXT TIME LEVEL 
PARRAY - ARRAY OF PRESSURES FOR PLOTTING 
PRI - INITIAL PRESSURE RATIO ACROSS THE SHOCK 
PLTCNT - COUNTER FOR GRAPHICS ROUTINES 
Q - ABSOLUTE FLUID VELOCITY 

QARRAY - ARRAY OF VELOCITIES FOR PLOTTING 
QLBD - INITIAL VELOCITY AT LEFT BOUNDARY 
Qll - INITIAL VELOCITY LEFT OF THE DIAPHRAGM 
QRBD - INITIAL VELOCITY AT RIGHT BOUNDARY 
QRI - INITIAL VELOCITY RIGHT OF THE DIAPHRAGM 
QQ - Q+A^^S (EXTENDED RIEMANN VARIABLE) 

RBDDR - RIGHT BOUNDARY DENSITY RATIO 
RBDDRI - INITIAL RIGHT BOUNDARY DENSITY RATIO 
RBDPR - RIGHT BOUNDARY PRESSURE RATIO 
RBDPRI - INITIAL RIGHT BOUNDARY PRESSURE RATIO 
RBDPRS - VALUE OF 0 DENOTES A CONSTANT PRESSURE AT RIGHT 
BOUNDARY, WHILE 1 IS FOR ADJUSTABLE PRESSURE 
RBDTR - RIGHT BOUNDARY TEMPERATURE RATIO 
RBDTRI - INITIAL RIGHT BOUNDARY TEMPERATURE RATIO 
RBNDRY - DENOTES RIGHT BOUNDARY OPEN OR CLOSED 
RNODE - ARRAY FOR RIGHT MOST NODE TO BE CORRECTED IN CORRCT 
RR - Q-A^S (EXTENDED RIEMANN VARIABLE) 

S - ENTROPY 

SARRAY - ARRAY OF ENTROPY FOR PLOTTING 
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C SIGMA - SPATIAL LOCATION OF OISCCNTINUITIES 
C SIGMA(L,J) WHERE L INCICATES THE TYPE OF 

C DISCONTIWITY AND J INDICATES THE 

C TIME LEVEL! 1 - CURRENT LEVEL 

C 2 - LEVEL BEING CALCULATED 

C SK - INTEGER THAT DENOTES RELATIVE LOCATION OF SHOCK NEAR 
C BOUr40ARIES 

C SKIP - VARIABLE WHICH INDICATES HOW MANY TIME STEPS BETWEEN 
C CALLS TO OUTPUT ROUTINES 

C T - TIME SINCE INITIAL CONDITIONS 

C TRI - INITIAL TEMP. RATIO ACROSS THE SHOCK 

C VHEAD - VELOCITY OF THE HEAD OF THE EXPANSION 

C WAVE FOR THE EXACT SOLUTION. 

C VTAIL - VELOCITY OF THE TAIL OF THE EXPANSION 
C WAVE FOR THE EXACT SOLUTION. 

C VCDE - VELOCTIY OF THE CONTACT DISCONTINUITY 
C FOR THE EXACT SOLUTION. 

C VS - THE SHOCK SPEED( POSITIVE RIGHT, NEGATIVE LEFT) 

C VSE - VELOCITY OF THE SHOCK FOR THE EXACT 

C SOLUTION. 

C W - MACH NO. RELATIVE TO A STANDING SHOCK 

C XARRAY - ARRAY OF SPATIAL POSITIONS FOR PLOTTING 

C XEXACT - ARRAY OF SIX X VALUES FOR THE EXACT SOLUTION. 

C XINIT - INITIAL POSITION OF DISCONTINUITY FOR 

C EXACT SOLUTION PLOTTING. 

C X2 - LOCATION OF NODE TO RIGHT OF DISCONT. 

C ALONG THE SPACIAL AXIS. 

C Y - (N + D/2 

C YEXACT - ARRAY OF SIX DENSITY VALUES FOR THE EXACT SOLUTION. 

C 

C OTHER VARIABLES ARE DEFINED IN THE 

C SUBROUTINES WHERE THEY ARE USED 

C 

c ♦ ♦ 

C ♦ PROBLEM SET - UP ♦ 

C ♦ ♦ 

C ♦♦♦♦♦♦ 

C 

C THE PARTICULAR PROBLEM FOR THIS VERSION IS: 

C 

C SHOCK TUBE, SINGLE CENTERED DIAPHRAGM WITH 

C HIGH PRESSURE SIDE TO THE RIGHT. 

C 

C BOUNDARY CONDITIONS - LEFT END CLOSED, RIGHT END OPEN 

C 

C VARIABLE DECLARATIONS 

C 

DIMENSION I2 ( ^ ) ,X2( 4 ) ,XEXACT( 6 ) ,YEXACT( 6 ) ,LNODE ( 4 ) ,RNODE( 4 ) , 
C SIGMA(4,2) 



USER INPUT REQUIRED HERE 
- SET THE DIMENSIONS EQUAL TO N 



DIMENSION A( 101 ) ,Q( 101 1 ,QQ( 101 ) ,RR( 101 ) ,S( 101 ) , 

C NEWRRt 101 ) ,NEWS( 101 ) ,NEWQQ( 101 ) , 

C PARRAY( 101 ) ,DARRAY( 101 ) ,SARRAY( 101 ) , 

C QARRAY( 101 ) , XARRAY ( 101 ) 

INTEGER I , J ,N , JSTOP ,Y .GRAPHS .COUNT . PLTCNT .BDRY ,SK . RBNDRY . 

C SKIP .12 .LWPRES.HALT.LNODE .RHODE .LBNDRY .LBDPRS.RBDPRS 

DOUBLE PRECISION TRI .PRI .QLI .QRI .ORI .G .G1 .G2 .SIGMA.EE .NEWQG . 

C DELT.H.ND.Xa.AR.H.DQ.VS.T.A.Q.QQ.RR.S.NEWRR.NEWS. 

C LBOPRI .LBDTRI .LBDDRI .LBDPR .LBDDR .LBDTR .GLBD . 

C RRA.QQA.AA.SA.QA.RRB.QQB.AB.SB.GB. 

C RBDPRI .RBDTRI .RBDDRI .RBDPR .RBDDR .RBDTR.GRBD . 

C QCS.VTEW.VHEW 

REAL VTAIL.VCDE.VSE.DLSH.DLCD.XINIT, VHEAD, 
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C XEXACT , YEXACT , PARRAY , OARRAY ,SARRAY >QARRAY ,XARRAY 

COMMON AR,DQ,VS,H 

ENTER THE APPROPRIATE VALUES BELOH 

N= 101 
GRAPHS=1 

FOR GRAPH = 1 OR a MUST ENTER CHANGES IN SUBROUTINE 

BORDER AND SUBROUTINE EXACT 

LINES ‘'FIRST ORDER N = ???** 

"OENSITY RATIO = ??? TEMP RATIO = ???'* 

“PRESSURE RATIO = ???“ 



SKIP=18 

JSTOP=101 

TRI=1.DDD0 

PRI=5.DD00 

DRI=5.DDD0 

QLI^D.DDDD 

QRI=0.0DDD 

LBDPRI = l.DDO 

LBDTRI = l.DOO 

LBDDRI = l.DDD 

RBDDRI = l.DOO 

RBDPRI = 4.DDD/5.D00 

RBDTRI = l.DOO 

G=1.4DD0 

EE=0.1D-a 

DENOTE LOW PRESSURE SIDE BY SETTING 

LWPRES = 2 IF LOW PRESSURE ON RIGHT 

LWPRES = 3 IF LOW PRESSURE ON LEFT 

LWPRES = 3 

SET BOUNDARY CONDITIONS BY SPECIFING OPEN OR CLOSED 

LBNORY: OPEN = 0, CLOSED = 1 ( FOR LEFT BOUNDARY) 

IF OPEN SPECIFY IF PRESSURE IS TO BE MAINTAINED AT LBDPRI 

OR IF IT CAN ADJUST TO PREVENT ANY WAVES FORMING AT THE BOUNDARY - 

LBDPRS: CONSTANT = 0, ADJUSTABLE = 1 

DO the same for RIGHT BOUNDARY > RBNDRY ,RB0PRS 

LBNDRY = 1 
LBDPRS = 1 
RBNDRY = 0 
RBDPRS - D 

exact SOLUTION VALUES 

XINIT=D.5D 

VHEAD= l.D 

VTAIL= 0.310557 

VCDE=-. 57^487 

VSE=-1. ^023^6 

DLCD=2. 713115 

DLSH=1. 69344 

SIGMA( 1»2)=D.5DDDDDD1D00 

SIGMA( 2»2)=D.5DDODDD1DOD 

SIGMAt 3,2 )=D .5DDDDD01DDD 

SIGMA(4,2)=D.5DDDDDD1DD0 

♦+++♦♦♦♦♦♦ END OF USER INPUT AREA 

SK = D 
T=0.0DD0 
DO ID 1=1,4 
12(1 )=0 
LNODEd ) = 0 
RNODE(I) = 0 
X2m=D.0De0 
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SIGMA(I,1)=0.0D00 

10 CONTirWE 

ND-DBLE(N) 

H=1.000/(ND-1.000) 

DO 11 1=1, N 

A(1 )=0.0D00 
Q(I)=0.DD00 
NENQQd )=0.0D00 
NENRRd )=0.0D0D 
NEWSd)=0.0D00 
QARRAYd)=D.D 
PARRAYd )=D.O 
0ARRAY(I)=0.0 
SARRAYC I )=0.0 

XARRAY( I ) = FL0AT( I-l ) mSNGL( H ) 

11 CONTINUE 

DELT=2.0000 

— LOAD INITIAL REIMAN VALUES INTO NODE LOCATIONS, FIRST FROM NODE— 

1 TO MIDPOINT (Y), AND THEN FROM Y TO N. NOTE IF SHOCK DOES NOT- 

—START AT MIDPOINT THEN Y SHOULD BE SET TO NODE WHERE SHOCK 

—INITIALLY IS 

AR=1.0D00 

DQ=D.0D00 

W=1.DD00 

VS=0.0D00 

QCS=0 . DOO 

VHEW=0,D00 

VTEW=0.D00 

G1=1.D00/(G-1.D00) 

G2=2,D00/(G-1.D00 ) 

— FLOW RIGHT 

IF(LWPRES,£Q.2) THEN 
LBDDR = DRI * LBDDRI 
LBDPR = PRI * LBDPRI 
LBDTR = TRI ^ LBDTRI 
QLBD = QLI 
RBDDR = RBDDRI 
RBDPR = RBDPRI 
RBDTR = RBDTRI 
QRBD = QRI 

Y = (N+D/2 
Y = 38 

DO 12 1=1, Y 

S( I )=G2-( 61/G )i^DLOG( PRI/( ( DRI ) ) 

QQC I )=QLI+DSQRT( TRI )»St I ) 

RR( I )=QLI-0SQRT( TRI )^^S( I ) 

12 CONTINUE 

Y = (N+3)/2 

Y=39 

DO 13 I=Y,N 
Sd )=G2 

QQ(I)=QRI+Sd) 

RRdl=QRI-Sd) 

13 CONTINUE 
ELSE 



FLOW LEFT 



LBDDR 


= 


LBDDRI 






LBDPR 


= 


LBDPRI 






LBDTR 


= 


LBDTRI 






QLBD ■ 


r 1 


QLI 






RBDDR 


= 


RBDDRI 




DRI 


RBDPR 


= 


RBDPRI 




PRI 


RBDTR 


= 


RBDTRI 




TRI 
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c 



QRBO s QRI 

Y » (N-D/2 
Y * 13 

DO I=1,Y 
Sm=G2 
QQm=QLI+S(I) 

RR(I)=QLI-S(I) 

14 CONTINUE 

Y = (N+D/2 

C Y = 14 

DO 15 I=Y,N 

S( I )=G2-( Gl/G )i(DLOG( PRI/( ( DRI ) ) 

QQi I )=QRI+DSQRT( TRI )^Si I ) 

RR( I )=QRI-DSQRT(TRI )^Si I ) 

15 CONTINUE 
END IF 

SET UP GRAPHICS PLOTS OF VARIABLES 

IF ( GRAPHS. GT.O) THEN 
CALL COMPRS 
CALL TEK618 

CALL HWROn 'AUTO' ) 

CALL HHSCALl 'SCREEN* ) 

IF ( GRAPHS. EQ. 2) THEN 
CALL PAGE( 11.0,8.5) 

ELSE 

CALL PAGE(8.5,11.0) 

END IF 

IF (GRAPHS. EQ.l) THEN 
CALL BORDER( JSTOP ) 

END IF 
END IF 
HALT = 0 
J=1 

COUNT =1 

IF (GRAPHS. EQ.l) THEN 

CALL PLOT(J, JSTOP, N,QQ,RR,S,H,XARRAY,PARRAY, 
ftOARRAY ,QARRAY ,SARRAY ,G ,G1 ,G2 ) 

END IF 

BURST DIAPHRAGM 

CALL TIME(N,QQ,RR,S,OELT,H) 

CALL 0BURST(N,H,QQ,RR,S,G,G1,G2,DELT,I2,X2,W,AR,DQ,VS,LWPRES, 

C SIGMA, A, Q) 

IF ( GRAPHS. EQ.O) THEN 

CALL LIST( N, SIGMA, QQ,RR,S,G,G1,G2,J,T,DELT, VS, QCS, VIEW, VHEH, 
C XINIT , VHEAD , VTAIL , VCOE ,VSE ) 

END IF 



BEGIN CALCULATION FOR JUMP TO NEXT TIME AND CONTINUE 

UNTIL EITHER JSTOP REACHED OR SHOCK MEETS CONTACT SURFACE 

16 IF (J.EQ. JSTOP) GOTO 18 
PLTCNT=J/SKIP 
C 

CALL TIME(N,QQ,RR,S,DELT,H) 

C 

CALL TRAK(N,SIGMA,H,QQ,RR,S,G,G1,G2,DELT,I2,X2,W,AR,0Q,VS,J, 

C LWPRES , QCS , VHEH ,VTEW ) 

C 

CALL SWEEP(N,H,SIGMA,QQ,RR,S,DELT,EE,Q,A,NEWQQ,NEWRR,NEWS,I2,G2, 
C J , LNOOE , RNODE , HALT , LBNDRY , LBOPRS , LBOPR , LBOTR , LBOOR , 

C QLBO ,G ,G1 ,RBNDRY ,RBDPRS ,RBDPR ,RBDTR ,RBODR ,QRBD ,3DRY , 

C SK) 

C 

IF(LNOOE(4).LT.J) THEN 
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U40DE(1) « 0 
RNODEdi s 0 
GO TO 17 

END IF 

IF(RN0DE(<f).LT.LN0DEC4n THEN 
RNODE(l) = 0 

END IF 
C 

CALL CORRCT(LNODE>RNODE>N>SIGHA,H>QQ»RR»S»G,GI,G2,iafX2>W>AR>DQ> 

C VS,A,Q) 

C 

17 IF (SIGMA(1,2KGE.1.D00) THEN 
IF(RBN0RY.EQ.01 THEN 

IF(SIGMA(1,2KNE.3.D00) THEN 
SK = 2 

CALL BONDRY(Q(N),Q(N-ll,QRBD,A(N),A(N-l),QQ(N),QQ(N-l), 
C RR(N),RR(N-1),S(N),S(N-1),H,EE,DELT, 

C RBNDRY ,RBDPRS >RBOPR ,RBOOR ,RBDTR , J ,NEWQQ( N ) , 

C NEWRR( N ) ,NEWS( N ) ,G ,G1 ,G2 >HALT >BDRY >SK ) 

END IF 

ELSE 

CALL EXTRAP(RR(N-11,RR(N-2),QQ(N-1),QQ(N-2),S(N-1), 

C S( N-2 ) >H >H ,RRA ,QQA ,SA ,AA »QA ) 

CALL SRFLCT( QQA ,RRA ,SA , SIGMA , VS ,OELT , LHPRES , 

C RR(N)>QQ(N)>S(N)»Q(N),A(N),G»G1»G2) 

END IF 

END IF 
C 

IF (SIGMAa,2).LE.O.OOO) THEN 
IF( LBNORY.EQ.O ) THEN 

IF(SIGMA( 1,2).NE.-2.D00) THEN 

BDRY = 3 

SK=2 

CALL BONDRYi Q( 1) 2 ) ,QLBD » A( 1) , A( 2 1 >QQ( 1 ) fQQi 2 ) , 

C RR( 1),RR( 2),S( 1 ),S( 2),H,EE,0ELT, 

C LBNDRY,LBDPRS,LBOPR>LBODR>LBDTR»J,NEHQQ( 1)» 

C NEHRR( 1 ) ,NEHS( 1 ) ,G,G1 ,G2,HALT ,BDRY ,SK ) 

END IF 

ELSE 

CALL EXTRAPt RR( 2 ) >RR( 3 ) ,qQ( 2 ) »QQ( 3 ) ,S( 2 ) > 

C S(3)>H>H,RRB,QQB,SB»AB,QB ) 

CALL SRFLCT(QQB,RRB,SB,SIGMA,VS,DELT,LWPRES, 

C RRl 1)»QQ( 1 ),S( 1 J »Qi 1 ),A( 1 ),G,G1>G2) 

END IF 

END IF 

T=T-fDELT 

C 

C OITTPUT DATA 

C 

IF ( ( COUNT . EQ . PLTCNTJ^SKIP ) . AND . ( GRAPHS . EQ . 1 ) ) 

C THEN 

C IF( (J.GT.SS)) THEN 

CALL PLOT( J , JSTOP ,N ,QQ ,RR ,S , H ,XARRAY ,PARRAY , 

CDARRAY ,QARRAY ,SARRAY ,G,G1 ,G2 ) 

C END IF 

END IF 

IF ( ( COUNT , EQ . PLTCNT*SKIP ) . AND . ( GRAPHS . EQ . 0 ) ) 

C THEN 

CALL LIST( N, SIGMA, QQ,RR,S,G,G1,G2,J,T,DELT, VS, QCS,VTEH,VHEH, 

C XINIT,VHEAD,VTAIL,VCOE,VSE ) 

END IF 

I F H COUNT . EQ . P LTCNT^SKI P ) . AND . ( GRAPHS . EQ . 2 ) ) 

C THEN 

IF (J.GT.50) THEN 

CALL EXACT! N,XINIT ,T ,VHEAO ,VTAIL ,VCDE ,VSE ,DLCO ,DLSH ,QQ ,RR ,S,H , 
CXARRAY ,0 ARRAY ,G ,G1 ,G2 ,ORI 1 

END IF 

END IF 
C 



149 



s 



o o o o o o o o o o o o o oooo 



TF(HAa.EQ.X) GO TO Id 
J=J+l 

CaJNT=COUNT-H 
GO TO 16 
Id CALL DONEPL 
END 

SUBROUTINE LIST( N ,SIGMA ,QQ ,RR ,S ,G ,G1 ,G2 , J ,T ,DELT ,VS ,QCS,VTEH, 
C VHEW,XINIT,VHEAD,VTAIL,VCDE,VSE ) 

♦ •fr + 4 + ♦♦ + ♦♦ + + ♦ + ♦ + ♦ + ♦♦ + ♦♦♦♦♦♦ 

+ 4 

♦ TABULAR RESULTS SUBROUTINE + 

♦ + 



VARIABLE DEFINITIONS 



DENS - DENSITY 
PRESS - PRESSURE 
TEMP - TEMPERATURE 

INTEGER I,J,N,L 

DIMENSION SIGMA( 4,2 ) ,QQ( N ) ,RR( N 1 ,S( N ) 

DOUBLE PRECISION SIGMA, QQ,RR,S, PRESS, VTEH,VHEW,QCS,TENXE, 
C TEMP, DENS, G,Gl,G2,Q,r,0ELT, VS, HEWXE , 

C XINIT ,VHEAD , VTAIL , VCOE , VSE ,3KXE ,CSXE 



TIME LEVEL*, J,' ELAPSED TIME IS*,T 

TIME STEP IS',0ELT,* SHOCK VELOCITY IS', VS 
CONTACT SURFACE VELOCITY ISSQCS 
HEAD EXPANSION WAVE VELOCITY IS',VHEW 
TAIL EXPANSION WAVE VELOCITY IS*,VTEW 



NODE 



VELOCITY 



DENSITY 



WRITE( 9,J^) 

WRITE( 9,^^) 

WRITEC 9,3^) 

WRITE( 9, 3*) 

WRITE! 9,3t) 

WRITE! 9, 3^) 

WRITE! 9, ») 

CPRESSURE* 

WRITE!9,3^) * * 

DO 61 1 = 1, N 

TEMP=( QQ( I )-RR( I ) )^i QQ( I )-RR! I ) )/( 4,DOO*S! I )*Si I ) ) 

DENS=( ! 1, 000/TEMP )3tDEXP! G3^! l.DOO-G)»(S! I )-G2)))»3M-Gl) 

PRESS=TEMP3^0ENS 

Q=(QQ!I)+RR(I))/2,ODOO 

WRITE (9,65) I, Q, DENS, PRESS 

65 FORMAT ! 4X,I2,7X,F12.6 ,7X,F12.6,7X,F12,6 ) 

61 CONTINUE 

WRITE(9,3t) • * 

WRITE! 9, • • 

WRITE! 9, 3^) • NODE QQ 

CMODIFIEO S* 

WRITE! 9, * * 

DO 62 1=1, N 

WRITE (9,66) I,QQ!I),RR!I),S!I) 

66 FORMAT !4X,I2,7X,F12.6,7X,F12,6,7X,F12.6 ) 

62 CONTINUE 

WRITE! 9, ») 

WRITE!9,3t) 

WRITE! 9,3^) 

WRITE! 9,3^) 

WRITE! 9,* ) 

DO 63 L=l,4 

WRITE!9,^0 L,* 

63 CONTINUE 

WRITE! 9,3^) • • 

WRITE!9,») * 

C 

WRITE! 9, * * 



RR 



DISCONTINUITY LOCATIONS AT TIME LEVEL *,J 



TYPE 



LOCATION* 



* , SIGMA! L, 2) 
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WRITE(9,M) • • 

IF(J.EQ.l) THEN 

WRITE ( 10, * TIME SHOCK CONTACT HEAD 

C TAIL* 

WRITE t 10, * • 

END IF 

WRITE (10,67) T ,$IGMA( 1,1) ,SIGMA( 2,1) ,SIGMA( 3 , 1 ) ,SIGMA( 4,1) 

67 FORMAT ( 1X,F12.6,1X,F12,6,1X,F12.6,1X,F12.6,1X,F12.6 ) 

IF(J.EQ,1) THEN 

WRITE(8,)() * EXACT VALUES* 

WRITE ( 8, ^^) • TIME SHOCK CONTACT HEAD 

C TAIL* 

WRITE(8,)() • * 

END IF 

SKXE=XINITWSE*T 

csxe=xinitwcdej^t 

HEWXE=XINITWHEAD^T 
TEWXE sXINITWTAIL^^T 
WRITE (8,68 ) T,SKXE ,CSXE ,HEWXE , TEWXE 

68 FORMAT ( 1X,F12. 6 ,1X,F12. 6 ,1X,F12. 6 , 1X,F12. 6 ,1X,F12. 6 ) 

RETURN 

END 

SUBROUTINE TIME( N ,(?t3,RR,S,DELT ,H ) 






CALCULATE TIME STEP SUBROUTINE 






NEW VARIABLE DEFINITIONS 

TMIN - RUNNING VALUE OF THE MINIMUM TIME STEP 



INTEGER N,I 

DIMENSION QQ(N),RR(N),S(N) 

DOUBLE PRECISION H,A,QQ,RR,S,DELT,TMIN,Q 

TMIN=2.0D00 

DO 21 1 = 1, N 

A=(QQ( I )-RR( I ) )/( 2.D00*S( I ) ) 
Q=(QQ(I)+RR(I))/2.0D00 
DELT=H/( DABS( DABS( Q )+A ) ) 

IF (DELT.LT.TMIN) THEN 
TMIN^OELT 
END IF 
21 CONTINUE 

DELT=0.99D00»TMIN 

RETURN 

END 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



SUBROUTINE TRAK(N,SIGMA,H,QQ,RR,S,G,G1,G2,DELT,I2,X2,W,AR,DQ,VS, 
«J,LWPRES,QCS,VHEW,VTEW) 



* DISCONTINUITY TRACKING SUBROUTINE * 



VARIABLE DE FINITIONS 



CSDIR - CONTACT SURFACE DIRECTION, 2 TO THE RIGHT, 3 TO THE LEFT 
CSRMN"- RIEMANN VARIABLE CHANGE ACROSS A CONTACT SURFACE 
DR - THE RATIO OF THE DENSITY ACROSS A 
SHOCK, B/A( RIGHT ),3/A( LEFT) 

DREMN - DUMMY VARIABLE 

PR - THE RATIO OF THE PRESSURE ACROSS A 
SHOCK, A/B( RIGHT ),B/A( LEFT) 



151 



ooooono oooo ooo oooonooooono 



MREim • THE MEASURED JUMP IN Qq ACROSS THE SHOCK, 

FROM A TO B, 

EREIMN - THE JUMP IN QQ ACROSS THE SHOCK CALCULATED ANALYTICALLY 
AS A FUNCTION OF H. 

SAP - ENTROPY TO THE LEFT OF THE SHOCK FOR FLOWS RIGHT 

OR ENTROPY TO THE RIGHT OF THE C.S. FOR FLOWS LEFT 
SBP - ENTROPY TO THE RIGHT OF THE C.S. FOR FLOWS RIGHT 

OR ENTROPY TO THE LEFT OF THE SHOCK FOR FLOWS LEFT 
SHKDIR - SHOCK DIRECTION OF TRAVEL, 3 TO THE LEFT,C TO THE RIGHT 
TS - TIME FOR SHOCK TO TRAVEL ONE INTERVAL 
X - DISTANCE FROM LEFT BOUNDARY TO NODE 

INTEGER N,I,Y,I2,L,J,SHKDIR,CSDIR,LWPRES 
DIMENSION SIGMA(4,2),X2(4),RR(N),QQ1N),$(N),I2(4) 

DOUBLE PRECISION SIGMA ,X2 ,X ,H ,AB ,SA , SB ,AA ,QA ,QB,QQA ,QQB ,RRA ,RRB , 
C RR,QQ,S,TS,W,0Q,AR,PR,G,G1,G2,VS,0ELT,CSRMN, 

C Q, MREIMN,DREMN, EREIMN, WW,SA1,SA2, SAP, SBP 

C AW ,QW , VHEW , VTEW ,TIME ,QCS 

LOCATING THE UPSTREAM NODE 

DO 10 L=l,4 

SI GMA ( L , 1) =SI GMA ( L , 2 ) 

Y=0 

X=0.D00 

1=1 

11 IF ( •NOT,(Y.EQ.On GOTO 10 

IF (SIGMAa,l).LT.X) THEN 
X2(U=X 
I2(L) = I 
Y=1 

END IF 
X=X+H 
1=1 + 1 
GOTO 11 
10 CONTINUE 

IF SHOCK HAS LEFT AN OPEN BOUNDARY OUT OF THE TUBE THEN SET 

SHOCK TO NEUTRAL 

IF (I2( D.GT.N) THEN 

SIGMA(1,1) = 2.DDD 
SIGMA(1,2) = 3.DDD 
W = l.DDO 



VS 


= 1.000 


DQ 


= 0.000 


AR 


s 0.000 


PR 


= l.DOO 


DR 


= 1.000 


GO 


TO 150 



ELSE IF(I2(1).LT.2) THEN 
SIGMA! 1,1) =-l.DD0 
SIGMA! 1,2) =-2.DDD 
W = l.DDD 



VS 




•l.DDD 


DQ 


r 


D.DDD 


AR 


r 


D.DDD 


PR 


= 


l.DDD 


DR 


= 


l.DDD 


GO 


TC 


1 15D 



END IF 

♦++++♦♦♦♦++++ DETERMINING SHOCK SPEED +♦++++♦♦ 

IF! ! J.EQ.D.OR. !I2! l).HQ.2).OR.!I2! D.EQ.N)) THEN 

-AT TIME ZERO OR BOUNDARYS DETERMINE CORRECT SHOCK DIRECTION 
-SHKDIR = 3 IS A SHOCK HEADED LEFT, AND SHKDIR = 2 IS SHOCK- 
-TRAVELING RIGHT 
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c 



IF(LWPRES.EQ,5) THEN 
SHKDIR = 3 
IF (J.EQ.l) THEN 

X2(l) = X2C1) - H 

12 ( 1 ) = 12 ( 1 ) - 1 

X2(2) = X2(2) - H 

12 ( 2 ) = 12 ( 2 ) - 1 

END IF 
GO TO 20 

ELSE 

SHKDIR = 2 

END IF 
GO TO 20 

—IF SHOCK AND CONTACT SURFACE ARE NOT WITHIN THE SAME 

INTERVAL THEN NO CORRECTIONS ARE NEEDED IN CALCULATING 

—THE REIMAN VARIABLE JUMP ACROSS THE SHOCK 

ELSE IF (12(1). NE. 12(2)) THEN 

IF ( ( SHKDIR . EQ . 3 ) . AND . ( SIGMA( 1 , 1 ) . EQ . ( X2 ( 1 )-H ) ) ) THEN 
X2(l) = X2(l) - H 
12 ( 1 ) = 12 ( 1 ) - 1 

END IF 
GO TO 20 

—IF SHOCK AND CONTACT SURFACE ARE WITHIN THE SAME INTERVAL 

THEN CORRECTIONS ARE REQUIRED TO DETERMINE SHOCK STRENGTH 

ELSE IF(SHKDIR.EQ.3) THEN 

—SHOCK LOCATION RELATIVE TO THE CONTACT SURFACE FOR A SHOCK 

—HEADED TO THE LEFT DETERMINES THE CORRECT VALUES FOR W AND VS— 

IF(SIGMA( 1,1).GT.SIGMA( 2>1) ) THEN 
GO TO 111 

ELSE IF(SIGMA(1,1).EQ.(X2(1)-H)) THEN 
X2( 1) = X2( 1) - H 
12 ( 1 ) = 12 ( 1 ) - 1 
X2(2) 3 X2(2) - H 
12 ( 2 ) = 12 ( 2 ) - 1 

END IF 

IS RRA=RR(I2(1)-1) 

RRB=RR(I2(D) 

QQA=QQ(I2(1)-1) 

QQB=QQ(I2(1) J 
SA=S(I2(1)-1) 

QA =(QQA+RRA)/2.D00 
QB =(QQB+RRB )/2.D00 
AA =(QQA-RRA)/( 2.D00)^SA) 

DQ =(QB-QA)/AA 

W = DSQRT( (DQJ^^*2)*(0.36D00) + 1.D00) - (DQ*0.6D00) 
DQ =(-l.D00)i(DQ 
GO TO 110 

FOR SHOCK HEADED RIGHT THE SAME PROCEDURE FOR CORRECTIONS ARE 

—FOLLOWED 

ELSE 

ELSE 
17 



IF(SIGMA(1,1).LT.SIGMA(2,1)) THEN 
GO TO 111 

RRA=RR(I2(1)-1) 

RRB=RR(I2(1) ) 
QQA=QQ(I2(1)-1) 
QQB=QQ(I2( 1) ) 

SB=S(I2( 1) ) 

QA =(QQA4RRA)/2.D00 
QB =(QQB4RRB )/2.D00 
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AB =(QQB-RRB)/(2,D00^^SB) 

DQ =(QA-QB)/AB 

W s DSQRTt (DQ)^*2)?^(0.36D00) + 1,D00) ♦ t0Q*0.6D00) 
GO TO 110 

ENO IF 

WITH NO SHOCX/CONTACT SURFACE INTERACTION THE JUMP IN REIMAN 

VARIABLES ARE DETERMINED WITHOUT INTERPOLATION OVER THE 

INTERVAL. MREIMN IS THE MEASURED JUMP, EREIMN IS THE ANALYTICAL— 

VALUE 

20 RRA=RR(I2(1)-1) 

RRB=RR(I2(D) 

QQA=QQ(I2(1)-1) 

QQB=QQ( 12(D) 

SA=S(I2(1)-1) 

SB-S( 12(D) 

21 AB=(QQB-RRB)/(2.D00^SB) 

AA=( QQA-RRA )/( 2 . DOO^SA ) 

IF(SHKDIR.EQ.3) THEN 

MREIMN = (RRB-RRA)ZAA 
DREMN = 0ABS( MREIMN) 

ELSE 

MREIMN = ( QQA-Q(5B )/AB 
DREMN = MREIMN 

END IF 

ITERATE FOR PROPER VALUE OF W USING THE QUADRATIC FIT OF THE 

REIMAN VARIABLE CHANGE WITH W CURVE. NOTE LEFT MOVING SHOCKS 

ARE USED IN THESE EQUATIONS SINCE RRB-RRA/AA=-( QQA-QQB/AB ) 

100 WW= (3.0396408001-((DREMN+2,7574D00)/0.286337D00)) 
W=5.51329A000-0SQRT(WW) 

DQ=2.000»(W^W-1.D00 )/( W^^( G+1 . DOO ) ) 

AR=0SQRT( 2.000*(G-1.000 1.000 + ( (G-l.DOO )*W^W/2 . 000 ) )^ 

C (G^G2*W)^W-1.000 ) )/( (G+1.000 J*W ) 

PR=( 2.D00*G/1G+1.000 ) ))^W*W-( ( G-1 . 000 I/(G+1.D00) ) 

0R=1 (G-1.000 I^WJ^W+2.000 )/( (G+1.000 )*W^^W ) 

ERElMN=0Q+( AR-1. DOO )i^G2-( AR^Gl/G )>^OLOG( PR>^( DR^^*G ) ) 

IF 1DABS(EREIMN-DABS( MREIMN )).LT. 0.10-5) GO TO 110 
OREMN = (OABSl MREIMN) - EREIMN) + OREMN 
GOTO 100 

SHOCK VELOCITY DEPENOS ON OIRECTION SHOCK IS TRAVELING 

LEFT IS < 0, ANO RIGHT IS > 0 

110 IF (J.EQ.D THEN 

TIME = O.DOO 

SAl = (Gl/G)M0LOG((2.D00J^GJ^(W^e^2)-G>1.000)/(G+l.D00)) 

SA2 = Gl)^0LOG(((G-1.000)J^(W*^2)+2)/((G+1.000)}^(W**2)n 
IF(SHK0IR.EQ.2) THEN 





SAP 


= SB - 


SAl - SA2 




SBP 


= SAP 




ELSE 


SBP 


= SA - 


SAl - SA2 




SAP 


= SBP 




ENO IF 









103 IF(SHK0IR.EQ.2) THEN 

CSRMN =( (DEXP( (SBP-SA)/G2) )i^(SA)-(SBP) )^^AR 

ELSE 

CSRMN =( (DEXP( (SAP-SB)/G2))i^(SB)-(SAP))*AR 

ENO IF 

IF(SHKDIR.EQ.3) THEN 

MREIMN =( (RRB-RRA)/AA)+CSRMN 
OREMN = DABS( MREIMN) 

ELSE 

MRE IMN =( ( QQA-QQB )/AB ) -CSRMN 
DREMN = MREIMN 

END IF 
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101 HW=( 3 . 0396403001-( ( 0REI1N+2 . 7574000 )/0 . 286337D00 ) ) 

H=5.513294DOO-OSORT( WW ) 

DQ=2 • OOOJf ( WJ^H-1 . DOO )/( mt G-M • DOO ) ) 

AR=DSQRT ( 2 , 000^( G-1 , 000 \*i 1 , 000+ ( ( G-1 • 000 , 000 ) )» 

C (G*G2*W*W-1.D00 ))/( (G+1,D00)*W) 

PR=( 2 . D00ifG/( G+1 . 000 ) )^tW»W-( ( G-1 . 000 )/( G+1 • 000 ) ) 

0R=( ( G-1 . DOO )JfW^^W+2 . 000 )/( ( G^l . DOO ))^W*W ) 

EREIMN=DQ+( AR-1 • 000 )^G2-( AR^^Gl/G )?^DLOG( PR^i DR^f^^G ) ) 

IF (DABS(EREIMN-DABS(MREIMN)).LT, 0.10-5) GO TO 102 ^ 
DREMN = (DABS(MREIMN) - EREIMN) ♦ DREMN 
GOTO 101 

102 SAl = <G1/G)*OLOG( (2.D00*G*(W**2)-G+1.D00)/(G+1.D00)) 

SA2 = G1^0L0G(MG-l.D00)*IN)^*2) + 2)/((G4l.D00)i^(W**2))) 
IF(SHK0IR.EQ.2) THEN 

SAP = SB - SAl - SA2 

ELSE 

SBP = SA - SAl - SA2 

END IF 

IF (OABS(SAP-SBP).LT. 0.10-5) GO TO 105 
IF(SHKDIR.EQ.2) THEN 

SBP = (SAP -SBP) + SBP 

ELSE 

SAP s (SBP-SAP) + SAP 

END IF 
GO TO 103 

END IF 

105 IF(SHK0IR.EQ.3) THEN 

DQ = ( -1.000 )i^D0 

VS = ( (RRA+QQA)*0.S000) - (W*AA) 

ELSE 

VS=( QQB+RRB )^^0 . SDOO+W^^AB 

END IF 

TS=H/DABS( VS) 

111 IF (TS.LT.DELT) THEN 
DELT=0.99000*TS 
END IF 

SIGMA( 1,2)=VS*DELT+SIGMA( 1,1) 

♦♦♦♦♦♦DETERMINE CONTACT SURFACE SPEE0++++^+^+^ 

IF(J.EQ.l) THEN 

CSOIR = SHKDIR 

END IF 

CONTACT SURFACE TRAVELING RIGHT 

150 IF(CSDIR.EQ.2) THEN 

CONTACT SURFACE MOVING RIGHT, CHECK FOR SHOCK IN INTERVAL 

and calculate speed OF CONTACT SURFACE AS APPROPRIATE 

IF(J.EQ.l) THEN 

QB =( QQB+RRB )/2. 000 
QA = OQ^^AB ♦ QB 

CALL SKJUMP( AB,QB,SB,AR,0Q,VS,G,G1,W,AA,QA,SA) 

Q = QA 
GO TO 200 

END IF 

IF(X2(2).EQ.X2(D) THEN 

IF(SIGMA(2,1).LE.SIGMA(1,D) THEN 

Q = (QQ(I2(2)-1) ♦ RR(I2(2)-D) / 2.000 

ELSE 

Q s(QQ(I2(2)) ♦ RR(I2(2))) / 2.000 

END IF 

ELSE 

QA = (QQ(I2(2)-1) ♦ RR(I2(2)-1)) / 2.000 
QB = (QQ(I2(2)) ♦ RR(I2(2))) / 2.000 
Q s (QA ♦ QB)/2.D00 

END IF 
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c 

c 

c 



c 

c 

c 



c 



c 

c 

c 

c 

c 

c 

c 

c 

c 



ELSE IF(CSDIR.EQ.3) THEN 

CONTACT SURFACE TRAVELING LEFT 

IF(J.EQ.l) THEN 

QA =(QQA+RRA)/2.DOO 
QB = DQ)(AA QA 

CALL SKJUMP ( AA ,QA >SA > AR , OQ , VS ,G ,G1 ,W , AB ,QB ,SB ) 

Q = QB 
GO TO 200 

END IF 

IF(X2(2).EQ.X2(D) THEN 

IF(SIGMA(2,1).GE.SIGHA( 1,1)) THEN 

Q = (QQ(I2(2n + RR(I2(2))) / 2.D00 

ELSE 

Q = (QQ(I2(2)-1) + RR(I2(2)-D) / 2.000 

END IF 

ELSE IF(SIGMA(2,1).EQ.<X2(2)-H)) THEN 
X2(2) = X2(2) - H 
12 ( 2 ) = 12 ( 2 ) - 1 
q = (QQ(I2(2))+RR(I2(2))) / 2.000 

ELSE 

QA = (QQ(I2(2)-1) ♦ RR(I2(2)-D) / 2.000 
QB = (QQ(I2(2)) ♦ RR(I2(2))) / 2.000 
Q = (QA ♦ QB)/2.D00 

END IF 

END IF 

200 SIGMA(2,2) = DELT M Q ♦ SIGMA(2,1) 

QCS=Q 

CALCULATE EXPANSION HAVE SPEED 

TIME = TIME+DELT 
IF(J.EQ.3) THEN 

AW=( QQ( ( N+1 )/2 )-RR( ( N+1 )/2 ) )/( 2 . DQO^Si ( N+1 )/2 ) ) 

QW=Q 

IF( (CSDIR).EQ.2) THEN 
VHEW=-( AW) 

VTEW= QW-AW 

ELSE 

VHEW= AW 
VTEW= QW+AW 

END IF 

END IF 

IF(J.EQ.3) THEN 

SIGMA(3,2) = SIGMA(3,1) + VHEW^TIME 

SIGMA(4,2) = SIGMA(4,1) ♦ VTEW^^TIME 

ELSE IF((CSDIR.EQ.2).AND.(SIGMA(3,1).GT.0.D00)) THEN 
SIGMA(3,2) = SIGMA(3,1) ♦ VHEW^DELT 

SIGMA(^,2) = SIGMA(4,1) ♦ VTEW*DELT 

ELSE IF( (CSDIR.EQ.3).AND.(SIGMA(3,1).GT.1.D00)) THEN 
SIGMA(3,2) = SIGMA(3,1) ♦ VHEW^DELT 

SIGMA(^,2) = SIGMA(4,1) ♦ VTEW^^DELT 

END IF 
RETURN 
END 

SUBROUTINE SWEEP(N, H, SIGMA, QQ,RR,S, DELT, EE ,Q, A, NEWQQ,NEWRR, NEWS, 
CI2,G2,J,LNODE,RNODE,HALT,LBNDRY,LBDPRS,LBDPR,LBDTR,LBDDR,QLBD,G, 
CG1,RBNDRY,RB0PRS,RB0PR,RB0TR,RB00R,QRBD,B0RY,SK) 

^ K 

it SPACE SWEEPING SUBROUTINE 

ik » 

)f**3«€*^****»^**9******9**^)t****^**************i(* 

VARIABLE DEFINITIONS 
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C AAVG - AVERAGE SPEED OF SOUND 

C CNTACT - I DIGIT VARIABLE DEWTING CONTACT SURFACE 
C LOCATICN>OIRECTION OF TRAVEL > AND IF IT CROSSES A NODE 

C DELQQH - CHANGE IN QQ FROM I TO I+l 

C DELQQL - CHANGE IN Oq FROM I-l TO I 

C DELRRH - CHANGE IN RR FROM I TO I+l 

C DELRRL - CHANGE IN RR FROM I-l TO I 

C DELSH - CHANGE IN S FROM I TO I+l 

C DELSL - CHANGE IN S FROM I-l TO I 

C DELAH - CHANGE IN A FROM I TO I+l 

C DELAL - CHANGE IN A FROM I-l TO I 

C DELQH - CHANGE IN Q FROM I TP I+l 

C DELQL - CHANGE IN Q FROM I-l TO I 

C DELX - INTERPOLATION DISTANCE ( LMD*DELT ) 

C DLTAit* - PREFIX WHICH INDICATES THE SPATIAL 

C CHANGE IN FOR ONE TIME STEP. 

C INTEG(K)- RESULT OF INTEGRATING Z(K) 

C JtHiNT - VALUE OF INTERPOLATED BETWEEN NODES 
C ON THE CURRENT TIME LEVEL. 

C LXX - NODE DEFINING THE LEFT INTERVAL 
C ^tPRIM(K)- SUFFIX WHICH INDICATES THE SPATIAL 
C DERIVITIVE OF * AT THE CURRENT TIME LEVEL. 

C RXX - NODE DEFINING THE RIGHT INTERVAL 
C SAVG - AVERAGE ENTROPY 

C SHOCK - 3 DIGIT VARIABLE DENOTING SHOCK 

C LOCATION, DIRECTION OF TRAVEL, AND IF IT CROSSES A NODE 

C JWSTEP - THE CHANGE IN TIME OF AT A NODE 
C USED TO STEP UP TO THE NEXT TIME LEVEL 

C X - LOCATION IN SPACIAL PLANE (I-l)J^H 

C Z(K1 - RIGHT SIDE OF THE K’TH EQUATION. 

C 

INTEGER I , RXX , LXX , SHOCK , CNTACT , 12 , J , LNODE , RNODE , HALT , 

C N » LBNDRY , LBOPRS ,RBNDRY ,RBDPRS ,SK ,BDRY 

DIMENSION SIGMAI4,2),S(N),Q(N),A(N),I2(<^),QINT(3),AINT(3),Z(3), 
C NEWQQ( N ) ,NEWRR( N ) ,NEWS( N ) .INTEG( 3 ) . APRIM( 3 ) ,QPRIM( 3 ) , 

C LNODE( ) ,RNODE( A ) ,AAVG( 3 ) ,RR( N ) ,QQC N ) 

DOUBLE PRECISION AAVG, SAVG, G2,X,H, SIGMA, QQ,RR,S,Q, A, G,G1, 

C DELQQH , DELQQL , DELRRH , DELRRL ,DELSH , DELSL , 

C DELAH, DELAL, DELQH, DELQL, DELT, 

C QINT , AINT ,QQINT ,RRINT , SINT , EE , 

C NEWQQ ,NEWRR ,NEWS , LBDPR , LBDTR ,LBDDR ,QLBD , 

C RRSTEP,SSTEP,INTEG,Z,DLTAQQ,QQSTEP, 

C DLTARR ,DLTAS ,APRIM ,QPRIM , 

C RBDPR .RBDTR ,RBDDR ,QRBD 

COMMON AR,DQ,VS,W • 

C 

c COMPUTE VELOCITY AND SPEED OF SOUND AT EACH NODE 

C 

DO 10 1= 1,N 

Q(I) = (QQ(I) ♦ RR(D) / 2.0D00 

Ad) = ((QQ(I) - RR(D) /(2.0D00 * S( I) ) ) 

10 CONTINUE 
C 

C ADVANCE LEFT BOUNDARY TO NEW TIME STEP 

C 

BDRY = 3 

IF(I2(1).EQ.2) THEN 
SK = 1 

ELSE IF(SIGMA(1,2).EQ.-2.D00) THEN 
SK = 3 

ELSE 

SK - 0 

END IF 

CALL BONDRY( Q( 1 ) ,Q( 2 ) ,QLBD , A( 1) , A( 2 ) ,QQ( 1 ) ,QQ( 2 ) ,RR(1 ) ,RR( 2 ) , 



C Sd ) ,S( 2 ) ,H ,EE , DELT , LBNDRY ,LBDPRS .LBDPR .LBDDR .LBDTR , 

C J ,NEWQQ( 1 ) ,NEWRR( 1 ) ,NEWSC 1 ) ,G ,G1 ,G2 .HALT .BDRY ,SK ) 

C 

C AT EACH NODE FROM 2 TO N-1 DETERMINE THE BEST ALGORITHM TO USE— 

C TO ADVANCE THAT NODE TO THE NEXT TIME STEP 
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c 

1=2 

11 IFd.EQ.Nl GO TO 1200 

X=FL0AT(I-1)J*H 
DELQQH = QQlI+l) - QQ(I) 

DELQQL = QQCI) - QQ( I-l) 

DELRRH = RR(I41) - RR(I) 

DELRRL = RR(I) - RR(I-l) 

DELSH s S(I+1) - S(I) 

DELSL = S(I) - S(I-l) 

DELAH = A(I + 1) - Ad) 

DELAL = ACI) - Ad-1) 

DELQH = Qd + 1) - Qd) 

DELQL = Qd) - Qd-1) 

DEFINE LEFT SECTOR AND RIGHT SECTOR HRT NODE EXAMINED— 

RXX = I 4 1 
LXX = I 

TEST FOR SHOCK 

IF (l2d).EQ.RXX) THEN 
SHOCK = 20D 
GO TO 20 

ELSE IF d2(l).EQ.LXX) THEN 
SHOCK = SOD 
GO TO 2D 

ELSE 

SHOCK = IDD 

END IF 
GO TO 3D 

DETERMINE DIRECTION SHOCK IS TRAVELING 

2D IF (SIGMAd,l).LT.SIGMA(l,2)) THEN 
SHOCK = SHOCK 4 20 

ELSE 

SHOCK = SHOCK + 3D 

END IF 

DETERMINE IF SHOCK CROSSES A NODE IN THIS TIME INTERVAL 

IF ( SHOCK. EQ. 220) THEN 

. IF (SIGMAC 1,2),GE.(X4H) ) THEN 
SHOCK = SHOCK 4 1 

ELSE 

SHOCK = SHOCK 4 2 

END IF 

ELSE IF ( SHOCK. EQ. 230) THEN 

IF (SIGMA( 1,2).LE.X) THEN 
SHOCK = SHOCK 4 1 

ELSE 

SHOCK = SHOCK 4 2 

END IF 

ELSE IF C SHOCK. EQ, 320) THEN 

IF (SIGMA(1,2).GE.X) THEN 
SHOCK = SHOCK 4 1 

ELSE 

SHOCK = SHOCK 4 2 

END IF 

ELSE IF ( SHOCK, EQ, 330) THEN 

IF (SIGMAd,2),LE,(X-H)) THEN 
SHOCK » SHOCK 4 1 

ELSE 

SHOCK = SHOCK 4 2 

END IF 

END IF 
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Q test for contact surface 

c 

30 IF <I2(2).EQ.RXX) THEN 
CNTACT s 200 
GO TO 40 

ELSE IF (I2(2).EQ.LXX) THEN 
CNTACT = 300 
GO TO 40 

ELSE 

CNTACT = 100 

END IF 
GO TO 50 
C 

C DETERMINE DIRECTION CONTACT SURFACE IS TRAVELING 

C 

40 IF (SIGMA(2,1).LT.SIGMA(2,2)) THEN 
CNTACT = CNTACT ♦ 20 

ELSE 

CNTACT 3 CNTACT ♦ 30 

END IF 
C 

C DETERMINE IF CONTACT SURFACE CROSSES A NODE DURING THIS TIME — 

C INTERVAL 

C 

IF C CNTACT. EQ. 220) THEN 

IF ISIGMA(2,2).GE.(X+H)) THEN 
CNTACT = CNTACT + 1 

ELSE 

CNTACT = CNTACT ♦ 2 

END IF 

ELSE IF ( CNTACT. EQ. 230) THEN 

IF (SIGMA(2,2).LE.X) THEN 
CNTACT = CNTACT ♦ 1 

ELSE 

CNTACT = CNTACT ♦ 2 

END IF 

ELSE IF (CNTACT.EQ.320) THEN 

IF (SIGMA(2,2).GE.X) THEN 
CNTACT = CNTACT ♦ 1 

ELSE 

CNTACT = CNTACT + 2 

END IF 

ELSE IF ( CNTACT. EQ. 330) THEN 

IF (SIGMA(2,2).LE.(X-H)) THEN 
CNTACT 3 CNTACT ♦ 1 

ELSE 

CNTACT = CNTACT + 2 

END IF 

END IF 
C 

C CHECK IF EITHER A SHOCK OR CONTACT SURFACE HITHIN H OF THIS NODE 

C DETERMINE PROPER ALGORITHM TO USE FOR CALCULATING EXTENDED REIMAN— 

C VARIABLE CHANGE ALONG CHARACTERISTICS AT THIS NODE 

C 

50 IF (SHOCK.EQ. 100. OR. CNTACT. EQ. 100) THEN 
C 

C NEITHER A SHOCK NOR A CONTACT SURFACE EXIST NEAR NODE— - 

C 

IF ( SHOCK.EQ. 100. AND. CNT ACT. EQ. 100) THEN 
C 

C TEST FOR SUBSONIC OR SUPERSONIC FLOH— 

C 

IF (0ABS(Q(I)).LT.A(D) THEN 

IF(I.EQ.(I2(l)-2)) THEN 
IF(J.EQ.l) THEN 

IF(SIGMA( 1,1).LE.SIGMA(1,2)) THEN 
GO TO 200 

END IF 

END IF 
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END IF 

IF(I.EQ.ri2(l)+l)) THEN 
IF(J.EQ.l) THEN 

IF( SIGMAC 1 ,1 ) .GE .SIGMA( 1,2)) THEN 
GO TO 300 

END IF 

END IF 

END IF 
GO TO 100 

ELSE 

IF (SIGMA(2,1).LE.SIGMA(2,2)) THEN 
GO TO 200 

ELSE 

GO TO 300 

END IF 

END IF 

END IF 

SHOCK OR CONTACT SURFACE ON LEFT, HEADED RIGHT, NO NODES CROSSED— 

IF (SHOCK. EQ. 322. OR. CNTACT.EQ. 322) THEN 
IF (DABS(Q(D).LT.A(I)) THEN 
GO TO 300 

ELSE 

GO TO 500 

END IF 

END IF 

SHOCK OR CONTACT SURFACE ON LEFT, HEADED LEFT, NO NODES CROSSED— 

IF ( SHOCK. £Q. 332. OR. CNT ACT. EQ. 332) THEN 
GO TO 300 

END IF 

SHOCK OR CONTACT SURFACE ON LEFT, HEADED RIGHT, NODE IS CROSSED 

IF (SHOCK. EQ. 321. OR. CNTACT.EQ. 321) THEN 
GO TO 400 

END IF 

SHOCK OR CONTACT SURFACE ON LEFT, HEADED LEFT, NODE IS CROSSED— 

IF ( SHOCK. EQ. 331. OR. CNTACT.EQ. 331) THEN 
GO TO 300 

END IF 

SHOCK OR CONTACT SURFACE ON RIGHT, HEADED RIGHT, NO NODES CROSSED— 

IF (SHOCK. EQ. 222. OR. CNTACT.EQ. 222) THEN 
GO TO 200 

END IF 

—SHOCK OR CONTACT SURFACE ON RIGHT, HEADED LEFT, NO NODES CROSSED— 

IF (SHOCK. EQ. 232. OR. CNTACT.EQ. 232) THEN 
IF (DABS(Q(D).LT.A(I)) THEN 
GO TO 200 

ELSE 

GO TO 500 

END IF 

END IF 

—SHOCK OR CONTACT SURFACE ON RIGHT, HEADED RIGHT, NODE CROSSED— 

IF (SHOCK. EQ. 221. OR. CNTACT.EQ. 221) THEN 
GO TO 200 

END IF 

—SHOCK OR CONTACT SURFACE ON RIGHT, HEADED LEFT, JUMPS NODE.— 
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C 



GO TO 400 
END IF 



BRANCH HERE IF SHOCK TO RIGHT OF CONTACT SURFACE 

—DETERMINE PROPER ALGORITHM TO USE FOR CALCULATING EXTENDED REIMAN 
—VARIABLE CHANGE ALONG CHARACTERISTICS AT THIS NODE 

IF (SIGMA(1,1).GT.SIGMA(2,D) THEN 

SHOCK ON LEFT, HEADED RIGHT, NO NODE JUMPED 

IF ( SHOCK. EQ. 322) THEN 

IF (CNTACT.EQ.322) THEN 

IF (DABS(Qd) ).LT.A(D) THEN 
GO TO 300 

ELSE 

GO TO 500 

END IF 

ELSE IF (CNTACT.EQ.332) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.331) THEN 
GO TO 800 

ELSE 

GO TO 900 

END IF 

END IF 

SHOCK ON LEFT, HEADED LEFT, NO NODE JUMPED— 

IF (SHOCK. EQ. 332) THEN 

IF (CNTACT.EQ.322) THEN 
GO TO 300 

ELSE IF (CNTACT.EQ.332) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.331) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.321) THEN 
GO TO 700 

ELSE 

GO TO 900 

END IF 

END IF 

SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE — 

IF (SHOCK. EQ. 321) THEN 

IF (CNTACT.EQ.322. OR. CNTACT.EQ.321) THEN 
GO TO 400 

ELSE IF (CNTACT.EQ.332. OR. CNTACT.EQ.331) THEN 
GO TO 800 

ELSE 

GO TO 900 

END IF 

END IF 

SHOCK ON LEFT, HEADED LEFT, JUMPS NODE — 

IF (SHOCK.EQ.331). THEN 

IF (CNTACT.EQ.322. OR. CNTACT.EQ.321) THEN 
GO TO 700 

ELSE IF ( CNTACT.EQ.332. OR. CNTACT.EQ.331) THEN 
GO TO 600 

ELSE 

GO TO 900 

END IF 

END IF 

SHOCK ON RIGHT, HEADED RIGHT, NO NODE JUMPED— 
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IF ( SHOCK. EQ. 222 I THEN 

IF (CNTACT.EQ.222) THEN 
GO TO 200 

ELSE IF (CNTACT.EQ.232.0R,CNTACT,EQ.231) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.321.0R.CNTACT.EQ.322) THEN 
GO TO 400 

ELSE IF (CNTACT.EQ.332.0R.CNTACT.EQ.33n THEN 
GO TO 600 

ELSE 

GO TO 900 

ENO IF 

END IF 

-—SHOCK ON RIGHT, HEAOEO LEFT, NO NOOE JUMPE0-— 

IF (SHOCK. EQ. 232) THEN 

IF (CNTACT.EQ.222) THEN 

IF (SIGMA(1,2).LT.SIGMA(2,2)) THEN 
GO TO 700 

ELSE 

GO TO 200 

ENO IF 

ELSE IF (CNTACT.EQ.231.0R.CNTACT.EQ.232) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.221) THEN 
GO TO 700 

ELSE IF (CNTACT.EQ.322) THEN 
GO TO 400 

ELSE IF (CNTACT.EQ.332.0R.CNTACT.EQ.331) THEN 
GO TO 600 

ELSE IF (SIGMAC 1,2).LT.SIGMA( 2,2) ) THEN 
GO TO 710 

ELSE 

GO TO 400 

ENO IF 

ENO IF 

—SHOCK ON RIGHT, HEAOEO RIGHT, JLIMPS NOOE 

IF ( SHOCK. EQ. 221) THEN 

IF (CNTACT.EQ.221. OR. CNTACT.EQ.222) THEN 
GO TO 200 

ELSE IF (CNTACT.EQ.231.0R.CNTACT.EQ.232) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ. 321. OR. CNTACT.EQ.322) THEN 
GO TO 400 

ELSE 

GO TO 800 

ENO IF 

END IF 

—SHOCK ON RIGHT, HEAOEO LEFT, JUMPS NOOE — 

IF (CNTACT.EQ. 221. OR. CNTACT.EQ. 222) THEN 
GO TO 700 

ELSE IF ( CNTACT.EQ. 231.0R.CNTACT.EQ. 232) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.322) THEN 

IF (SIGMA(1,2).LT.SIGMA(2,2)) THEN 
GO TO 710 

ELSE 

GO TO 400 

END IF 

ELSE IF (CNTACT.EQ. 321) THEN 
GO TO 710 

ELSE IF (OABS(Q(I)).LT.A(in THEN 
GO TO 600 
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ELSE 



GO TO 800 



END IF 

—BRANCH HERE IF SHOCK IS TO LEFT OF CONTACT SURFACE 

—DETERMINE PROPER ALGORITHM TO USE FOR CALCULATING EXTENDED REIMAN 
VARIABLE CHANGE ALONG CHARACTERISTICS AT THIS NODE 

ELSE IF (SIGMA(l,l).LT.SIGMA(2,in THEN 

—SHOCK ON LEFT, HEADED RIGHT, NO NODE CROSSED— 

IF ( SHOCK. EQ. 122) THEN 

IF( CNTACT . EQ . 222 . OR . CNTACT . EQ . 221 ) THEN 
GO TO bOO 

ELSE IF (CNTACT.EQ.332) THEN 

IF (SIGMA( 1,2 ).GT.SIGMA(2,2)) THEN 
GO TO 700 

ELSE 

GO TO 300 

END IF 

ELSE IF (CNTACT.EQ.331) THEN 
GO TO 700 

ELSE IF (CNTACT. EQ. 222. OR. CNTACT. EQ. 221) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.232) THEN 
GO TO <►00 

ELSE IF (SIGMA(1,2).GT.SIGMA(2,2)) THEN 
GO TO 710 

ELSE 

GO TO 400 

END IF 

END IF 

—SHOCK ON LEFT, HEADED LEFT, NO NODE CROSSED— 

IF ( SHOCK. EQ. 332) THEN 

IF (CNTACT.EQ.322) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.332) THEN 
GO TO 300 

ELSE IF (CNTACT. EQ. 321) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.222.0R.CNTACT.EQ.221) THEN 
GO TO 800 

ELSE IF (CNTACT. EQ. 231. OR. CNTACT. HQ. 232) THEN 
GO TO 400 

ELSE 

GO TO 900' 

ENO IF 

END IF 

—SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE— 

IF ( SHOCK. EQ. 321) THEN 

IF (CNTACT. EQ. 322. OR. CNTACT. EQ. 321) THEN 
GO TO 600 

ELSE IF (CNTACT, EQ. 332. OR. CNTACT.EQ.331) THEN 
GO TO 710 

ELSE IF (CNTACT.EQ. 222. OR. CNTACT. EQ. 221 ) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ. 231) THEN 
GO TO 710 

ELSE IF (SIGMA(1,2).GT.SIGMA(2,2)) THEN 
GO TO 710 

ELSE 

GO TO 400 

ENO IF 

END IF 
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SHOCK ON LEFT, HEADED LEFT, JUMPS NODE™ 

IF (SHOCK. EQ. 331) THEN 

IF (CNTACT.EQ.322.0R.CNTACT.EQ.321) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.331.0R.CNTACT.EQ.332) THEN 
GO TO 300 

ELSE IF (CNTACT.EQ.221.0R.CNTACT.EQ.222) THEN 
GO TO 800 

ELSE 

GO TO 400 

END IF 

END IF 

—SHOCK ON RIGHT, HEADED RIGHT, NO NODE CROSSED — 

IF ( SHOCK. EQ. 222) THEN 

IF (CNTACT.EQ.222.0R.CNTACT.EQ.221) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.231) THEN 
GO TO 710 

ELSE IF (CNTACT.EQ.232) THEN 

IF (SIGMA(1,2).GT.SIGMA(2,2)) THEN 
GO TO 700 

ELSE 

GO TO 200 

END IF 

ELSE 

GO TO 900 

END IF 

END IF 

—SHOCK ON RIGHT, HEADED LEFT, NO NODES CROSSED — 

IF (SHOCK. HQ. 232) THEN 

IF (CNTACT.EQ.222.0R.CNTACT.EQ.221) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.232) THEN 

IF (DABS(Q(D).LT.A(I)) THEN 
GO TO 200 

ELSE 

GO TO 500 

END IF 

ELSE 

GO TO 900 

END IF 

END IF 

—SHOCK ON RIGHT, HEADED RIGHT, JUMPS NODE — 

IF ( SHOCK. EQ. 221) THEN 

IF (CNTACT.EQ.221.0R.CNTACT.EQ.222) THEN 
GO TO 600 

ELSE IF (CNTACT.EQ.231) THEN 
GO TO 710 

ELSE IF (CNTACT.EQ.232) THEN 
GO TO 700 

ELSE 

GO TO 900 

END IF 

END IF 

—SHOCK ON RIGHT, HEADED LEFT, JUMPS NODE — 

IF (CNTACT.EQ.221.0R.CNTACT.HQ.222) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.231. OR. CNTACT.EQ.232) THEN 
GO TO 400 



% 
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ELSE 



GO TO 900 



END IF 

—SHOCK AND CONTACT SURFACE ARE AT THE SAME LOCATION AFTER TIME 

ZERO — ■ 

ELSE IF( ( SHOCK. EQ. 222 )• AND. CCNT ACT ,EQ. 222 ) ITHEN 
SHOCK = 321 
CNTACT = 321 
GO TO 400 

ELSE IF( ( SHOCK . EQ . 322 ) . AND . ( CNTACT . EQ . 322 ) )THEN 
IF(DABS(Q(in.LT.A(D) THEN 
GO TO 300 

ELSE 

GO TO 500 

END IF 

ELSE IF( (SHOCK. EQ. 332). AND. (CNTACT. EQ. 332)) THEN 
SHOCK = 231 
CNTACT = 231 
GO TO 400 

ELSE IF( (SHOCK. EQ. 232). AND. (CNTACT. EQ. 232) ) THEN 
IF(DABS(Q(I)).LT.A(I)) THEN 
GO TO 200 

ELSE 

GO TO 500 

END IF 

ELSE 

GO TO 720 

END IF 

CALL CONDITION SUBROUTINE WHICH CONTAINS ALGORITHM THAT IS— 

NUMERICALLY BEST SUITED FOR THE SITUATION AT NODE I — 

100 CALL CONDl(Q(I),Q(I+l),A(I),A(I+l),RR(I),QQ(I),S(I),DELQQL, 

C DELRRH,DELSH,DELSL,DELQH,DELAH,OELQL,DELAL,H,EE, 

C DELT ,QQINT ,RRINT ,S1NT ,QPRIM , APRIM , AINT ,Q( I-l ) , A( I-l ) ) 

GO TO 1000 

200 CALL C0ND2(Q(I),Q(I-1),A(I),A(I-1),RR(I ),QQ(I),S(I),DELQQL, 

c oelrrl,delsl,delql,delal,delt,h,ee,qqint,rrint,sint, 

C QPRIM, APRIM, AINT) 

GO TO 1000 

300 CALL C0ND3(Q(I),Q(I+1),A(I),A(I>1),RR(I),QQ(I),S(I),DELQQH, 

C OELRRH,DELSH,OELQH,DELAH,DELT,H,EE,QQINT,RRINT,SINT, 

C QPRIM, APRIM, AINT) 

GO TO 1000 

400 CALL CON04(I,SHOCK,CNTACT,J,LNOOE,R^JODE,QQSTEP,RRSTEP,SSTEP) 

GO TO 1100 

500 CALL C0ND5(Q(I),Q(I-1),Q(I41),A(I),A(I-1),A(I^1),RR(I ),QQ(I), 

C S(I ),DELQQL,DELQQH,DELRRL,DELRRH,DELSL,OELSH,OELQL, 

C DELQH,DELAL,DELAH,H, EE, DELT,QQINT,RRINT, SINT, AINT, 

C QPRIM, APRIM, SHOCK, CNTACT) 

GO TO 1000 

600 CALL C0ND6( SHOCK, CNTACT, HALT) 

QQSTEP-0.00 
RRSTEP=0.D0 
SSTEP=0.D0 
GO TO 1100 

700 CALL CON07(SHOCK, CNTACT, DELT, SIGMA, 12, 1, H, HALT, Q(D) 

QQSTEP=Q.OO 
RRSTEP=O.DO 
SSTEP=0\D0 
GO TO 1100 
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c 

710 CALL COt^D7N(SHOCK,CNTACT,DELf>SIGMA,I2,I,H,HALT,Q(I),X) 
QQSTEP=0.D0 
RRSTEP=0.D0 
SSTEP=0.D0 
GO TO 1100 
C 

720 CALL CON07S( SIGMA, HALT, SHOCK, CNT ACT) 

QQSTEP=0.D0 
RRSTEP=0.D0 
SSTEP=0.D0 
GO TO 1100 
C 

300 CALL C0N08( SHOCK, cm* ACT, HALT) 

QQSTEP-O.DO 
RRSTEP=0.00 
SSTEP=0.00 
GO TO 1100 
C 

900 PRINT » ,*AN IMPOSSIBLE SITUATION EXISTS - ERROR' 
QQSTEP=0.D0 
RRSTEP=0.D0 
SSTEP=0.D0 
HALT = 1 
GOTO 1100 

CALCULATE DLTA QQ, DLTA RR i DLTA S 

1000 DLTAQQ=QQINT-QQ( I ) 

DLTARR=RRINT-RR(I) 

OLTAS=SINT-S(I) 

CALCULATE Z(K)'S 



AAVGC 1 ) = ( AINT( 1 )+A( I ) )/2 . 0000 
AAVG( 3 ) = ( AINTI 3 )+A( I ) )/2 . 0000 
AAVG( 2 )=0.0000 
SAVG s (SINT+S(I) )/2.000 

Z( 1)=-C 1 . 0000/G2 )^»AAVG( 1 )»( SAVG-G2 QPRIM( 1 )-G2»APRIM( 1 ) ) 
2( 3 ) = ( 1 . 0D00/G2 )»AAVG( 3 SAVG-G2 )«( QPRIMC 3 )+G2^f APRIM( 3 ) ) 

2( 2)=0.0D00 

INTEGRATE THE 2(K)'S 

INTEGt 2)=0.0000 

INTEG(l)=Za)J(OELT 

INTEG(3)=Z(3)i<0ELT 

SOLVE THE EQUATION 

QQSTEP=OLTAQQ+INTEGI 1 ) 

RRSTEP=DLTARR4lNTEG( 3 ) 

SSTEP=OLTAS*INTEG( 2 ) 

STORE THE SOLUTION 

1100 NEWQQ( I )=QQ( I )+QQSTEP 

NEWRRt I )=RR( I )4RRSTEP 
NEWS(I)=S(I)+SSTEP 

GO TO NEXT NOOE 

1=1 + 1 
GOTO 11 
1200 CONTINUE 

BORY = 2 

, IF(I2(1).EQ.N) THEN 

SK = 1 
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ELSE IF(SIGMAfl,2).EQ.5,000) THEN 
SK = 3 

ELSE 

SK » 0 

END IF 

CALL BONORY(Q(N)»Q(N-I)>QRBO»A(N),A(N-1),QQ(N)»QQ(N-1)>RR(N)> 
C RR(N-1),S(N),S(N-1),H,EE,DELT,RBNDRY,RBDPRS,RBDPR,RBDDR, 
C RBDTR » J ,NEHQQ( N ) >NEHRR( N ) »NEHS( N ) ,G1 ,G2 ,HALT >BDRY ,SK ) 

UPDATE THE VARIABLES 



1=1 

1210 IF (I.EQ.N+1) GOTO 1220 
RR(I)=NEWRR(I) 
QQ(I)=NEWQQ(I) 
S(I)=NEWS(I ) 

I=l4l 
GOTO 1210 
1220 CONTINUE 
RETURN 
END 









CONDITION SUBROUTINES 1-8 









SUBROUTINE CONDK QI ,QIP1 ,AI ,AIP1 ,RR ,QQ ,S ,DELQQL ,DELRRH ,DELSH , 

C DELSL>0ELQH,0ELAH,0ELQL,0ELAL,H,£E,0ELT,QQINT, 

C RRINT,SINT,QPRIM>APRIM,AINT»QIM1,AIM1) 

ik * 

^ SUBROUTINE CONDITION 1 * 



USED WHEN NO CONTACT SURFACES NOR SHOCKS WITHIN H OF NODE 
AND FLOW IS SUBSONIC 

CALCULATES QQINT , RRINT tSINT ,QPRIM , APRIM 

USES FORWARD AND BACKWARD DIFFERENCE SCHEMES 



VARIABLE DEFINITIONS 

AI - A(I) 

AIMl - Ad-D* 

AIPl _ A(I+1) 

E(K) - ACTUAL ERROR IN CHARACTERISTIC SLOPE CALCULATION. 

LMD - SLOPE OF THE CHARACTERISTICS IQ+A,Q,Q-A) 

QI - Q(I) 

QIMl - Q(I-l) 

QIPl - Q(I+1) 

DIMENSION LMD( 3 ) ,DELX( 3 ) ,QINT( 3 ) , AINT( 3 ) ,E ( 3 ) ,QPRIM( 3 ) ,APRIM( 3 ) 
INTEGER K 

DOUBLE PRECISION QI , QIPl tAl , AIPl ,RR ,QQ »S , AIMl , QIMl , 

C DELQQL ,DELRRH ,DELSH ,DELSL ,DELQH ,DELAH ,DELQL , 

C DELAL ,H ,EE ,DELT , LMD ,DELX ,QINT , AINT ,E , 

C QQINT tRRINT ,SINT ,QPRIM , APRIM 

INITIAL ESTIMATE OF CHARACTERISTIC SLOPES 

LMD(l) = QIMl 4 AIMl 
LMD(2) = QI 
LMD(3) = QIPl - AIPl 

CALCULATE LINEARLY INTERPOLATED VALUES OF Q AND A 
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C 



10 K s I 

20 IF (K.LT.4) THEN 

DELXIK) = OELT « LMO(K) 

IF (LMO(K).LT. 0.0000) THEN 
QI)rr(K) = QI - 



ELSE 



AINT(K) = AI - 



QINT(K) s QI - 



(OELX(K) 




BELQH 


/ 


H) 


(OELX(K) 




OELAH 


/ 


H) 


(DELX(K) 




DELQL 


/ 


H) 


(OELX(K) 




OELAL 


/ 


H) 



END IF 
K = K ♦ 1 
GO TO 20 



END IF 



CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE 
FROM NEW INTERPOLATED VALUES 

Ed) s OABS(LMDd) - (QINT(l) ♦ AINT(l))) 

E(2) = DABS(LMD(2) - (OINT(2))) 

E(5) = 0ABS(LMD(3) - (QINT(3) - AINT(3))) 



LMD(l) a QINTd) ♦ AINTd) 

LM0(2) a QINT(2) 

LM0(5) a QINT(3) - AINT(3) 

COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TIL MET 



IF (CEd).GT.EE).0R.(E(2).GT.EE).0R.(E(3),GT.EE)) GO TO 10 



CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN 

VARIABLES AND MODIFIED ENTROPY AT POINT A 

QQINT a QQ - (DELXd) » (OELQQL / H)) 

IF (LM0( 2). LE. 0.0000) THEN 

SINT a S - (DELX(2) ^ (DELSH / H)) 

ELSE 

SINT a S - (0ELX(2) * (OELSL / H)) 

END IF 

RRINT a RR - (0ELX(3) * (OELRRH / H)) 

CALCULATE SPATIAL DERIVATIVES — 



C 



C 

C 

C 

C 

C 

C 

C 

C 

c 

c 



QPRIMd) a 
QPRIM(2) a 
QPRIM(3) a 
APRIM(l) a 
APRIM(2) a 
APRIM(3) a 
RETURN 
END 



(QI-QIMD/H 

O.DOO 

(QIPl-QI )/H 
(AI-AIMD/H 
O.DOO 

( AIPl-AI )/H 



SUBROLTTINE CON02( QI ,QIM1 ,AI ,AIM1 ,RR ,QQ,S , OELQQL ,DELRRL , OELSL , 
C DELQL,DELAL,DELT,H, EE, QQINT, RRINT, SINT, QPRIM, 

C APRIM,AINT) 



* * 
^ SUBROUTINE CONDITION 2 * 



BACKWARD DIFFERENCE ALGORITHM FOR CALCULATING — 

RRINT, QQINT, SINT, APRIM, QPRIM, AINT 

DIMENSION LMD( 3 ) ,0ELX( 3 ) ,QINT( 3 ) ,AINT( 3 ) ,E ( 3 ) ,QPRIM( 3 ) ,APRIM( 3 ) 
INTEGER K 

DOUBLE PRECISION QI ,QIM1,AI, AIM1,RR,QQ,S, OELQQL, OELRRL, OELSL, 

C OELQL,DELAL, OELT, H, EE, AINT ,QINT,LMD,OELX,E, 
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s 



c 



QQirrr ,rrint ,sint ,qprim, aprim 



initial estimate of characteristic slopes 

LMD(l) = QIMl ♦ AIMl 
LMD(2) = QI 
LMD(3) = QI - AI 

CALCULATE LINEARLY INTERPOLATED VALUES OF Q AND A 

10 K s 1 

20 IF (K.LT.^^) THEN 

DELX(K) = DELT ^ LMD(K) 

QINT(K) = QI - (DELX(K) * DELQL / H) 

AINT(K) = AI - (DELX(K) * DELAL / H) 

K = K ♦ 1 
GO TO 20 

END IF 

CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE 

from new INTERPOLATED VALUES 

Ed) = DABS(LMDd) - (QINT(l) ♦ AINTd))) 

E(2) = 0ABS(LMD(2) - QINK2)) 

E(3) = DABS(LMD13) - (QINT(3) - AINT( 3 ) ) ) 

LMDCl) = QINTd) ♦ AINTd) 

LMD(2) = QINT(2) 

LMD(3) = QINT(3) - AINT(3) 

COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TIL MET 

IF ((Ed).GT.EE).OR.(E(2).GT.EE).OR.(E(3),GT.EE)) GO TO 10 

CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN 

VARIABLES AND MODIFIED ENTROPY AT POINT A 

QQINT = QQ - (DELXd) * (DELQQL / H)) 

SINT = S - (DELX(2) ^ (DELSL / H)) 

PRINT = RR - (DELX(3) * (DELRRL / H)) 

CALCULATE SPATIAL DERIVATIVES — 



QPRIMd) = 
QPRIM(2) = 
QPRIM(3) = 
APRIMd) = 
APRIM(2) = 
APRIM(3) = 
RETURN 
END 



(QI-QIMl )/H 
O.DOO 

(QI-QIMD/H 
( AI-AIMD/H 
O.DOO 

(AI-AIMD/H 



SUBROUTINE 

C 

c 



C0ND3(QI,QIP1,AI,AIP1,RR,(3Q,S,DELQQH,DELRRH,0ELSH, 
DELQH ,DELAH ,DELT ,H ,EE , QQINT , PRINT , SINT ,QPRIM , 
APRIM, AINT) 



* ^ 
SUBROUTINE CONDITION 3 * 



FORWARD DIFFERENCE ALGORITHM FOR CALCULATING 

PRINT, QQINT, SINT, APRIM, QPRIM, AINT 

DIMENSION LMD( 3 ) ,DELX( 3 ) ,QINT( 3 ) ,AINT( 3 ) ,E( 3 ) ,QPRIM( 3 ) ,APRIM( 3 ) 
INTEGER K 

DOUBLE PRECISION QI ,QIP1 ,AI ,AIP1 ,RR ,QQ ,S , 

C DELQQH ,DELRRH ,DELSH ,DELQH ,DELAH , 



non ooooooo o noon non non 



C DELT,H,EE>AINT,QIMT,LMD,OELX,E, 

C QQINT,RRIMT,SINT>QPRIH,APRIM 

INITIAL ESTIMATE OF CHARACTERISTIC SLOPES 

LMD(l) s QI 4 AI 
LM0(2) = QI 
LM0(5) = QIPl - AIPl 

CALCULATE LINEARLY INTERPOLATED VALUES OF Q AND A 

10 K s 1 

20 IF (K.LT,‘^) THEN 

DELX(K) = DELT LMOCKl 

QINT(K) = QI - (OELX(K) » DELQH / HI 

AINT(K) = AI - (DELX(K) » OELAH / H) 

K = K 4 1 
GO TO 20 

END IF 

CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE 

FROM NEW INTERPOLATED VALUES 



Ed) s DABS(LMDd) - (QINTd) 4 AINTd))) 

E(2) = DABS(LMD(2) - QINT(2)I 

E(3) s DABS(LMD(3) - (QINT(3) - AINT(3))) 

LMDd) = QINTd) 4 AINTd) 

LMD(2) = QINT(2) 

LMD(3) = QINT(3) - AINT(3) 

COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TIL MET— 

IF((Ed),GT.EE).OR.(E(2).GT.EE).OR.(E(3),GT,EE)l GO TO 10 

CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN 

VARIABLES AND MODIFIED ENTROPY AT POINT A 

QQINT s QQ - (DELXd) ^ DELQQH / H) 

SINT = S - (DELX(2) » DELSH / H) 

RRINT = RR - (DELX13) DELRRH / H) 

— CALCULATE SPATIAL DERIVATIVES — 



C 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 



QPRIMd) = 
QPRIM(2) = 
QPRIM(3) = 
APRIMd) = 
APRIMC2) = 
APRIMC3) = 
RETURN 
END 



(QIPl-QI )/H 
O.DOO 

(QIPl-QI )/H 
(AIPl-AIl/H 
O.DOO 

( AIPl-AI )/H 



SUBROUTINE C0ND4( I , SHOCK ,CNTACT , J , LNODE ,RNODE ,QQSTEP ,RRSTEP , 
C SSTEP ) 



* SUBROUTINE CONDITION 4 ^ 



A SITUATION EXISTS AT NODE I THAT WILL BE CORRECTED IN — 
SUBROUTINE CORRCT. ^^NODE ARRAYS ARE GIVEN INFORMATION — 
ON NODE LOCATION AND SHOCK/CONTACT SURFACE PICTURE 

- VARIABLE DEFINITIONS 

CD^TRK - DENOTES HOW MANY NODES NEED TO BE CORRECTED 
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THIS TIME STEP 

DIMENSION LNOOE(^),RNOOE(^) 

INTEGER LNODE, RHODE, I, SHOCK, CNT ACT >J,C04TRK 
DOUBLE PRECISION QQSTEP ,RRSTEP ,SSTE? 

-ASSIGN FIRST NODE ENCOUNTERED DURING THE NEW TIME STEP TO 
-LNOOE 

IF(LN0DE(4).LT.J) THEN 
CO^^TRK = 1 
END IF 

IFCCD^TRK.EQ.ll THEN 
LNODE(l) = I 
LNODE12) = SHOCK 
LNODE(5) = CNTACT 
LNODE(^) = J 



-IF A SECOND NODE WITH CONDITION 4 IS ENCOUNTERED IN THE 
-SWEEP, THIS NODE IS ASSIGNED TO RHODE 



ELSE 



END IF 



RNODE(l) = I 
RNODE(2) = SHOCK 
RNODE(5) = CNTACT 
RNODE(^) 3 J 



-LNODE AND RHODE WILL BE "JUMPED** OVER DURING SWEEP THUS- 
-THEIR **STEP VALUES ARE SET TO 0 

QQSTEP * O.DOO 
RRSTEP 3 O.DOO 
SSTEP 3 O.DOO 



IF(C04TRK.EQ.l) THEN 

CD4TRK s CD4TRK ♦ 1 

ELSE 

CD4TRK = 1 

END IF 
RETURN 
END 

SUBROUTINE C0ND5( QI ,QIM1 ,QIP1 , AI ,AIM1 ,AIP1 ,RR ,QQ ,S,DELQQL , 



C 

C 

C 



DELQQH,DELRRL,DELRRH>DELSL,DELSH,DELQL>DELQH, 
DELAL,DELAH,H,EE,DELT,QQINT,RRINT>SINT,AINT, 
QPRIM,APRIM, SHOCK, CNTACT ) 









SUBROUTINE CONDITION 5 






w wwvrwwww wwww wwwwwwwww w www’wwwwwwwwwwwwwwwwwwww wwwww wwwwww 



FOR SUPERSONIC FLOW WITH A DISCONTINUITY ON ONE SIDE OF THE NODE- 
CALCULATES QQINT,RRINT ,SINT ,APRIM,QPRIM,AINT — 

DIMENSION LMD( 3 ) ,DELX( 3 1 ,QINT( 3 ) ,AINT( 3 ) ,E( 3 ) ,QPRIM( 3 ) ,APRIM( 3 1 
INTEGER K, SHOCK, CNTACT 

DOUBLE PRECISION QI , AI ,RR,QQ ,S ,QIM1 ,QIP1 ,AIM1 ,AIP1 , 

C EE,DELT,LMD,DELX,QINT,AINT,E,DELQL,DELQH,DELAL, 

C QQINT ,SINT ,RRINT ,QPRIM ,APRIM ,DELAH ,H , 

C DELQQL,DELQQH,DELRRL,DELRRH,DELSL,DELSH 

DISCONTINUITY ON LEFT, HEADED RIGHT, NOT CROSSING NODE — 

IF(( SHOCK. EQ. 322 ). OR. I CNTACT. EQ. 322)1 THEN 

INITIAL ESTIMATE OF CHARACTERISTIC SLOPES . 
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LMD(l) = QI ♦ AI 
LMDC2) s QI 
IM0(5) = QIPl - AIPl 

CALCULATE LINEARLY INTERPOLATED VALUES OF Q AND A 

10 K s 1 

20 IF (K.LT.4) THEN 

DELX(K) = DELT ^ LMD(K) 

QINTCKI = QI - CDELX(K) » DELQH / H) 

AINT(K) = AI - (DELX(K) » DELAH / H) 

K = K ♦ 1 
GO TO 20 

END IF 

CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE 

FROM NEW INTERPOLATED VALUES 

Ed) = DABS(LMDd) - (QINT(l) ♦ AINT(l))) 

E(2) # DABS(LMD(2) - QINT(2)) 

E(3) » DABS(LMD(3) - (QINT(3) - AINT(5))) 

LMD(l) = QINT(l) ♦ AINT(l) 

LMD(2) = QINT(2) 

LMDCS) = QINTC3) - AINT(3) 

COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TO MEET 

IF( (E(1).GT.EE).0R.(E(2).GT,EE).0R.(E(5).GT.EE n GO TO 10 

CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN 

VARIABLES AND MODIFIED ENTROPY AT POINT A 

QQINT = QQ - (DELXd) ^ DELQQH / H) 

SINT = S - (DELX(2) * DELSH / H) 

RRINT = RR - (DELX(3) DELRRH / H) 

— - CALCULATE SPATIAL DERIVATIVES 

QPRIM(l) = (QIPl-QD/H 
QPRIM(2) = O.DOO 
QPRIM(3) = (OIPl-QD/H 
APRIM(l) = (AIPl-AD/H 
APRIM(2) = 0.000 
APRIM(2) = (AIPl-AD/H 
END IF 

DISCONTINUITY ON RIGHT, HEADED LEFT, NOT CROSSING NODE — 

IF( (SHOCK. EQ. 232). OR. (CNTACT.EQ. 232)) THEN 

INITIAL ESTIMATE OF CHARACTERISTIC SLOPES 

LMD(l) = QIMl ♦ AIMl 
LMD(2) = QI 
LMD(3) = QI - AI 

CALCULATE LINEARLY INTERPOLATED VALUES OF Q AND A 

30 K = 1 

^0 IF (K.LT.4) THEN 

DELX(K) = DELT * LMD(K) 

QINT(K) = QI - (DELX(K) » DELQL / H) 

AINT(K) = AI - (OELX(K) ^ DELAL / H) 

K = K ♦ 1 
GO TO 40 

END IF 

CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE 
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FROM NEM INTERPOLATED VALUES 

Ed) = DABS(LMDCl) - (QINT(l) ♦ AINT(l))) 

E(2) s DABS(LMD(2) - QINT(2)) 

E(S) = DABSCLMD(5) - (QINT(3) - AINT(3))) 

LMDd) = QINTd) ♦ AINT(l) 

LMD(2) = QINT(2) 

LMD(S) = QINT(5) - AINT(3) 

COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TO MEET 

IF ((Ed.).GT.HE).OR.(E(2).GT.EE).OR.(E(3).GT.HE)) GO TO 30 

CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN 

VARIABLES AND MODIFIED ENTROPY AT POINT A 

OQINT = QQ - (DELX(l) C DELQQL / H)) 

SINT = S - (DELX(2) ( DELSL / H)) 

RRINT = RR - (DELX(3) * (DELRRL / H)) 

CALCULATE SPATIAL DERIVATIVES 



QPRIM( 1) 

QPRIMC2) 

QPRIMC3) 

APRIM(l) 

APRIMC 2) 

APRIM(3) 

END IF 

RETURN 

END 



(QI-QIMl )/H 
O.DOO 

(QI-QIMl )/H 
( AI-AIMl )/H 
O.DOO 

(AI-AIMD/H 



SUBROUTINE COND6( SHOCK ,CNT ACT , HALT ) 



» * 
» SUBROUTINE CONDITION 6 ^ 



INTEGER SHOCK ,CNT ACT, HALT 

PRINT * ,'THE SITUATION FOR CONDITION 6 WITH THE CONTACT SURFACE * 
PRINT * ,'TO THE RIGHT OF A SHOCK BOTH HEADED RIGHT OR THE 
PRINT , ‘CONTACT SURFACE TO THE LEFT OF A SHOCK BOTH HEADED LEFT* 
PRINT » , ‘ADDITIONAL LOGIC IS REQUIRED TO PROCEED* 

PRINT ^ , 'SHOCK =*, SHOCK,* CONTACT SURFACE =',CNTACT 

HALT = 1 

RETURN 

END 



SUBROLTTINE COND7( SHOCK ,CNTACT ,DELT , SIGMA ,12,1 ,H ,HALT ,Q ) 

* SUBROUTINE .CONDITION 7 » 



DIMENSION SIGMA(^^,2),I2C4) 

INTEGER HALT, SHOCK, CNT ACT, 12, I 
DOUBLE PRECISION DELT , SIGMA ,H ,Q 

PRINT ^ , 'CONDITION 7 REQUIRES THAT THE CONTACT SURFACE AND ' 

PRINT » , 'SHOCK! MOVING IN OPPOSITE DIRECTIONS) MEET AND CROSS. * 

PRINT ^ ,*WHEN THEY INTERSECT , THE RESULT IS A FUNCTION OF EXIST-* 
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C 



PRINT » , 
PRINT ^ , 
PRINT ^ , 
PRINT ^ , 
PRINT ^ , 
PRINT , 
PRINT ^ , 
PRINT » , 
PRINT ^ , 
PRINT » , 
PRINT ^ , 
PRINT » , 
PRINT » , 
PRINT , 
HALT = 1 
RETURN 
END 



*ING CONDITIONS AROUND THEM. BOTH THE ORIGINAL SHOCK 
•AND CONTACT SURFACE WILL EXPERIENCE VELOCITY CHANGES . 
•THIS SUBROUTINE WOULD HAVE TO CALCULATE WHEN AND WHERE 
•WITHIN THE TIME/SPACE ITERVAL THE SHOCK AND CONTACT 
•SURFACE INTERSECT. THIS IS BASED ON KNOWN SPEED AND 
•THEIR RESPECTIVE LOCATIONS. THIS NEW TIME, OELT(NEH) 
•THEN COULD BE USED TO RERUN ''SWEEP^^, WITH CONDITION 7S 
'EXISTING AT TIME = T ♦ DELT(NEW). ADDITIONAL LOGIC IS 
•REQUIRED TO CONTINUE. • 

■SHOCK =•, SHOCK, • CONTACT SURFACE =^,CNTACT 
•SIGMA(1,1) =' ,SIGMA( 1,1), • SIGMA(1,2) =' ,SIGMA( 1,2 ) 
•SIGMA(2,1) =' ,SIGMA(2,l)f* SIGMA(2,2) = ' ,SIGMA( 2 ,2 ) 
•CONTACT SURFACE VELOCITY =* ,Q 
•NODE I =• ,I 



SUBROUTINE COND7N( SHOCK ,CNTACT ,DELT , SIGMA ,12 ,I ,H ,HALT ,Q,X ) 



^ ik 

^ SUBROUTINE CONDITION 7N 

ik ik 



DIMENSION SIGMA(4,2),I2(4) 

INTEGER HALT,SHOCK,CNTACT,I2,I 
DOUBLE PRECISION DELT ,SIGMA ,H ,Q,X 



PRINT « , 
PRINT » , 
PRINT » , 
PRINT » , 
PRINT * , 
PRINT ^ , 
PRINT » , 
PRINT ^ , 
PRINT * , 
PRINT * , 
PRINT » , 
PRINT * , 
PRINT » , 
PRINT ^ , 
PRINT * , 
PRINT * , 
PRINT ^ , 
PRINT , 
PRINT * , 
PRINT » , 
HALT = 1 
RETURN 
END 



•CONDITION 7N REQUIRES THAT WHEN EITHER A SHOCK OR 
•CONTACT SURFACE JUMPS A NODE (AS DETERMINED BY 
•COMPARING SIGMAtL,!) TO SIGMAIL,2)) A CONTACT SURFACE' 
•OR SHOCK WOULD BE MET AND CROSSED DURING THE JUMP. 

•THE RESULT WHEN THEY INTERSECT IS A FUNCTION OF 
•EXISTING CONDITIONS AROUND THEM. BOTH THE ORIGINAL • 
•SHOCK AND THE CONTACT SURFACE WILL EXPERIENCE VELOCITY' 
•CHANGES. THIS SUBROUTINE WOULD HAVE TO CALCULATE WHEN^ 
•WITHIN THE TIME INTERVAL AND WHERE SPACIALLY THE 
•INTERSECTION OCCURS. THIS IS BASED ON KNOWN SPEED AND* 
•THE RESPECTIVE LOCATIONS OF THE SHOCK AND CONTACT 
•SURFACE. THE NEW TIME, OELT(NEW), COULD THEN BE USED * 
•TO RERUN "SWEEP'^ WITH CONDITION 4 AT THIS NODE, AND 
•SIGMA(2,2) = SIGMA(1,2) SO THAT CONDITION 7S WOULD 
•RESULT IN THE NEXT TIME INTERVAL.' 

•SHOCK =• ,SHOCK, 'CONTACT SURFACE =',CNTACT 
•SIGMA(1,1) s*,SIGMA(l,l),* SIGMA(1,2) = * ,SIGMA( 1 ,2 ) 
•SIGMA(2,1) =',SIGMA(2,1),* SIGMA(2,2) = ' ,SIGMA( 2 , 2 ) 
•CONTACT SURFACE VELOCITY s' ,Q 
•NODE I =•,!,’ LOCATION AT X =*,X 



SUBROUTINE C0ND7S( SIGMA , HALT , SHOCK ,CNTACT ) 

^ * 

SUBROUTINE CONDITION 7S * 

ik » 



DIMENSION SIGMA(4,2) 

INTEGER SHOCK, CNTACT, HALT 
DOUBLE PRECISION SIGMA 

PRINT * , 'CONDITION 7S HAS BEEN MET. THIS MEANS THAT THE SHOCK ' 
PRINT » ,'AND CONTACT SURFACE ARE LOCATED AT THE SAME X AT A ' 
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PRINT 
PRINT 
PRItif 
PRirrr 
PRINT 
PRINT 
HALT a 1 
RETURN 
END 



.•TIME OTHER THAN ZERO. THIS SUBROITTINE HOULO HAVE TO 
.•DETERMINE THE RESULT OF THE INTERSECTION BASED ON THE 
.•CONDITIONS TO THE LEFT AND RIGHT. ADDITIONAL LOGIC IS 
.•REQUIRED TO PROCEED. • 

.•SHOCK =•, SHOCK.* CONTACT SURFACE a^.CNTACT 
.•SIGMAJl.l) a^.SIGMAU.l),’ SIGMA(2.1) =^,SIGMA(2.1) 



SUBROUTINE C0ND8( SHOCK. CNTACT, HALT) 

^ SUBROUTINE CONDITION 3 * 

INTEGER HALT .SHOCK, CNTACT 

PRINT ^ . 'CONDITION 3 COULD RESULT ONLY AFTER THE ORIGINAL SHOCK 
PRINT * ,'HAS CROSSED THE CONTACT SURFACE. AND THE SUBSEQUENT 
PRINT * . 'CONDITIONS DETERMINED. ADDITIONAL LOGIC IS REQUIRED TO 
PRINT » . 'CONTINUE.' 

PRINT .'SHOCK =•. SHOCK.* CONTACT SURFACE =*, CNTACT 

HALT a 1 

RETURN 

END 

SUBROUTINE CORRCT( LNODE .RNODE ,N .SIGMA ,H ,QQ .RR .S ,G tGl ,G2 ,12 ,X2 , 
C AR.DQ.VS.A.Q) 

« ^ 

^ DISCONTINUITY CORRECTION SUBROUTINE * 



VARIABLE DEFINITIONS 



ENTRPY - MODIFIED ENTROPY 
SNDSPO - SONIC VELOCITY 
UA - VELOCITY RELATIVE TO THE SHOCK, 
UB - VELOCITY RELATIVE TO THE SHOCK, 
VLCTY - VELOCITY 



LEFT SIDE 
RIGHT SIDE 



INTEGER LNODE .RNODE ,N ,12 .NODE .SHOCK .CNTACT ,K 

DIMENSION LNODE( 4 ) .RNODE ( 4 ) ,SIGMA( 4, 2 ) ,QQ( N ) ,RR( N ) .$( N ) ,I2( 4 ) . 
C X2(4),A(N),Q(N),XA(4).XBI4) 

DOUBLE PRECISION SIGMA ,QQ ,RR ,S ,A ,Q , 

C H,G,G1,G2,X2,W,AR,0Q.VS. 

C RR A , R RB . QQ A . QQB , A A , AB . Q A . QB , S A . SB . 

C SAl ,SA2 .VLCTY .SNDSPO . ENTRPY .QQCALC .RRCALC , 

C XB.XA 



-DEFINE STATEMENT FUNCTIONS- 



QQCALC( VLCTY .SNDSPO .ENTRPY ) 
RRCALCi VLCTY .SNDSPD .ENTRPY ) 



VLCTY + SNDSPD*ENTRPY 
VLCTY - SNDSPD*ENTRPY 



-SET INITIAL VALUES OF VARIABLES TO 2ER0- 



RRA a O.DOO 
RRB a 0.000 
QQA a O.DOO 
QQB a O.DOO 
AA a 0.000 
AB a O.DOO 
QA a 0.000 
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QB s O.DOO 
SA s O.DOO 
SB = O.DOO 
00 15 

XAtK)=0.00 

XB(K)=0.00 

15 CONTINUE 

DETERMINE IF ONLY ONE NODE NEEDS TO BE CORRECTED OR TWO 

NODES. ALSO, IF SHOCK AND CONTACT SURFACES ARE CLOSE 

ENOUGHtWITH IN 2H) TO INTERACT 

IF( ( LNODE ( 2 ) . EQ . 100 ) . OR . ( LNODE ( 3 ) . EQ . 100 ) ) THEN 

IF(( LNODE( 1 ) . GE . ( RNODE ( 1 )-2 ) ) . AND . ( RNODE ( 1 ) . GT . 0 ) ) THEN 
GO TO 20 

ELSE 

GO TO 30 

END IF 

ELSE IF(RNODE(1).GT.O) THEN 
GO TO 20 

ELSE 

GO TO 10 

END IF 

BRANCH HERE IF ONLY ONE NODE TO BE CORRECTED (LNODE) AND, 

SHOCK/CONTACT SURFACE INTERACTION 

10 SHOCK = LNODE(2) 

CNTACT = LNOOE(3) 

SHOCK ON LEFT, HEADED RIGHT, JUMPS OR NOTvCONTACT SURFACE ON 

—RIGHT, HEADED LEFT, DOESN’T CROSS NODE — 

IF( (SHOCK. EQ. 322. OR. SHOCK. EQ. 321). AND. (CNTACT. EQ. 232)) THEN 

DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE 

12(1) = LNODE(l) ♦ 1 
12(2) = LNODE(l) + 1 
CALL DELTAX( 12, SIGMA, XB,XA,H) 

E>aRAPOLATE TO RIGHT FACE OF CONTACT SURFACE 

IF (12(2). EQ.N) THEN 

CALL 3BDRY(RR(N),QQ(N),S(N),RRB,QQB,SB,AB,QB) 

ELSE 

CALL EXTRAP(RR(I2( 2)),RR(I2( 2)4l),QQ(I2( 2) ),QQ(I2( 2 )>1), 

C S(I2( 2)),S(I2( 2)4l),XB( 2),H,RRB,QQB,SB,AB,QB) 

END IF 

CALCULATE VARIABLE CHANGE ACROSS CONTACT SURFACE 

SA = S(I2(2)-1) 

CALL CSJUMP( AB,QB,SB,SA,G2,QA,AA) 

IF SHOCK INTERACTS CALCULATE CHANGE IN VARIABLES ACROSS SHOCK 

IF( SHOCK. EQ. 321) THEN 
QB = QA 
AB = AA 
SB = SA 

CALL SKJUMP( AB,QB,SB,AR,DQ,VS,G,G1,W,AA,QA,SA) 

QQA = QQCALC(QA,AA,SA) 

RRA = RRCALC(QA,AA,SA) 

INTERPOLATE TO NODE ( I2( 1 )-l ) AND ASSIGN CORRECTED VALUES 

IF (12(1). EQ. 2) THEN 

CALL BBDRY(RRA,QQA,SA,RR(I2(1)-1),QQ(I2(1)-1 ), 
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C 



C 

C 

C 



ELSE 



END IF 



S(ia(l)-1),A(I2(X)-1),Q(I2(1)-1)) 

ELSE 

CALL INTERPt RRA,RR( I2C 1 )-2 ) ,QQA,QQ( I2( 1 )-2 ) ,SA, 
SCI2( 1)-2),XA(1),(XA(1)+H),RR(I2C1)-1), 
QQ( I2C 1 )-l ) ,SC I2( 1 )-l ) ,A( I2C 1 )-l ) , 
Q(I2(1)-D) 

END IF 



QQ(I2f2)-l) - QQCALC(QA,AA>SA) 
RRCI2(2)-1) = RRCALC(QA,AA,SA) 
S(I2(2)-1) = SA 
A(I2(2)-1) = AA 
Q(I2(2)-1) = QA 



—SHOCK ON RIGHT, HEADED LEFT, CROSSES NODE OR NOTi CONTACT SURFACE- 
— ON LEFT, HEADED RIGHT, DOES'NT CROSS NODE 

ELSE IF( (SHOCK. EQ. 232. OR. SHOCK. EQ. 231). AND. CONTACT. EQ. 322)) THEN 

DETERMINE NODE TO RIGHT TO SHOCK AND CONTACT SURFACE 



12(1) = LNODE(l) 

12(2) = LNODE(l) 

CALL DELTAX(I2,SI(;MA,XB,XA,H) 

•EXTRAPOLATE TO LEFT FACE OF CONTACT SURFACE 

IF (12(2). EQ. 2) THEN 

CALL BBDRYC RR( 1 ) ,QQ( 1 ) ,S( 1 ) ,RRA ,QQA ,SA , AA ,QA ) 

ELSE 

CALL EXTRAP(RR(I2( 2)-l ),RR(I2( 2)-2),QQ(I2( 2)-l),QQ(I2( 2)-2), 
C S(I2( 2)-l),S(I2( 2)-2),XA(2),H,RRA,QQA,SA,AA,QA) 

END IF 

•CALCULATE VARIABLE CHANGE ACROSS CONTACT SURFACE 



SA = S(I2( 2) ) 

CALL CSJUMPC AA,QA,SA,SB,G2,QB,AB) 

•IF SHOCK INTERACTS CALCULATE CHANGE IN VARIABLES ACROSS SHOCK 

IF( SHOCK. EQ. 231) THEN 
QA = QB 
AA = AB 
SA - SB 

CALL SKJUMP( AA,QA,SA,AR,DQ,VS,G,G1,H,AB,QB,SB) 

QQB = QQCALC(QB,AB,SB) 

RRB - RRCALC(QB,AB,SB) 

•INTERPOLATE TO NODE(I2(D) AND ASSIGN CORRECTED VALUES 

IF (12(1). EQ.N) THEN 

CALL BBDRY( RRB ,QQB ,SB,RR( I2( 1 ) ) ,QQ( I2( 1 ) ) ,S( I2( 1 ) ) , 
A(I2(1)),Q(I2( 1))) 

ELSE 

CALL INTERP( RRB,RR( I2( 1 )♦! ) ,QQB ,QQ( I2( 1 ) + l ) ,SB, 

S(I2(1) + 1),XB(1),(XB(1)+H),RR(I2(1)), 

QQ(I2( 1)),S(I2(1)),A(I2(1)),Q(I2( 1))) 

END IF 

ELSE 

QQ(I2(2)) QQCALC(QB,AB,SB) 

RR(I2(2)) 3 RRCALC(QB,AB,SB) 

S(I2(2)) = SB 
A(I2(2)) = AB 
Q(I2(2)) s QB 

END IF 

- SHOCK ON RIGHT, HEADED RIGHT, DOESN'T CROSS OR ON LEFT, HEADED— 



C 

c 

c 
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RIGHT, DOES CROSS! AND CONTACT ON LEFT, HEADED RIGHT, CROSSES 

OR NOT: OR SHOCK ON RIGHT, HEADED LEFT, DOESN'T CROSS NODE 

HITH CONTACT ON LEFT, HEADED RIGHT, AND CROSSED NODE 

ELSE IF( ( (SHOCK. EQ. 222. OR. SHOCK. EQ. 321). AND. (CNT ACT. EQ. 321. OR. 

C CNT ACT. EQ. 322)). OR. ( SHOCK. E(9. 232. AND. CNTACT.EQ. 321)) 

C THEN 

DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE 

12(1) = LNODE(l) ♦ 1 
12(2) = LNODE(l) ♦ 1 
CALL DELTAX( 12, SIGMA ,XB,XA,H) 

EXTRAPOLATE TO RIGHT FACE OF SHOCK 

IF (12(1). EQ.N) THEN 

CALL BBDRY(RR(N),QQ(N),S(N),RRB,Q(3B,SB,AB,(^) 

ELSE 

CALL EXTRAP(RR(I2( 1)),RR(I2(1)+X),QQ(I2( 1) ),QQ(I2(1 )♦!), 

C S(I2(1)),S(I2(1)>1),XB(1),H,RRB,QQB,SB,AB,QB) 

END IF 

CALCULATE CHANGE IN VARIABLES ACROSS SHOCK 

IF( SHOCK. EQ. 232) THEN 
AA = AB/AR 

SAl = (Gl/G)»DLOG((2.DOO*G»(W^f^*2)-G>l.DOO) / (G+l.DOO)) 
SA2 = Gl)^DLOG( ( (G-l.DOO )>^(W**2)42.D00) / ((G^l.DOO)* 

C ( W^^*2 ) ) ) 

SA * SB>SA1+SA2 
UB = QB - VS 
UA = UB -AA»DQ 
QA = UA > VS 

ELSE 

CALL SKJUMP( AB ,QB ,SB ,AR ,DQ ,VS ,G ,G1 ,W , AA ,QA ,SA ) 

END IF 

CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE 

IF(CNTACT.EQ.322) THEN 

QQ(I2(1)-1) = QQCALC(QA,AA,SA) 

RR(I2(1)-1) = RRCALC(QA,AA,SA) 

S(I2(1)-1) = SA 
A(I2(1)-1) = AA 
Q(I2(1)-1) = QA 

ELSE 

QB = QA 
AB = AA 
SB = SA 

SA = S(I2(2)-2) 

CALL CSJUMP( AB ,QB ,SB ,SA ,G2 ,QA , AA ) 

QQA s QQCALC(QA,AA,SA) 

RRA s RRCALC(QA,AA,SA) 

INTERPOLATE TO NODE( I2( 2 )-l ) AND ASSIGN CORRECTED VALUES 

IF (12(2). EQ. 2) THEN 

CALL BBDRY( RRA ,QQA ,SA ,RR( I2( 2 )-l ) ,QQ( I2( 2 )-l ), 

C S(I2(2)-1),A(I2(2)-1),Q(I2(2)-D) 

ELSE 

CALL INTERP(RRA,RR( I2( 2 )-2 ) ,QQA ,QQ( I2( 2 )-2 ) ,SA , 

C S(I2( 2)-2),XA(2),(XA( 2)+H),RR(I2(2)-l), 

C QQ(I2( 2)-l),S(I2( 2)-l),A(I2( 2)-l), 

C Q(I2(2)-D) 

END IF 

END IF 

—SHOCK ON LEFT, HEADED LEFT, DOESN'T CROSS NODE, OR ON RIGHT,— 
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C —HEADED LEFT, CROSSED NODEj AND COmACT ON RIGHT, HEADED LEFT,— 

C —CROSSES NODE OR NOT: OR SHOCK 0T4 LEFT, HEADED RIGHT, DOESN'T— 

C CROSS NODE, AND CONTACT ON RIGHT, HEADED LEFT, CROSSED NODE — 

C 

ELSE 2F( ((SHOCK. EQ. 332. OR.SHOCK, EQ. 231). AND. (CNT ACT. EQ. 232. OR- 
C CNTACT . EQ . 231 ) ) . OR . ( SHOCK . EQ . 3 2 2 . AND . CKTACT . EQ . 231 ) ) 

C THEN 

DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE 

12(1) = LNODE(l) 

12(2) = LNODE(l) 

CALL DELTAX( 12, SIGMA, XB,XA,H) 

-EXTRAPOLATE TO LEFT FACE OF SHOCK 

IF (12(1). EQ. 2) THEN 

CALL BBDRY( RR( 1 ) ,QQ( 1 ) ,S( 1 ) ,RRA ,QQA ,SA ,AA ,QA ) 

ELSE 

CALL EXTRAP(RR(I2(l)-l),RR(I2tl)-2),QQ(I2(l)-l),QQ(I2(l)-2), 

C S(I2(1 )-l),S(I2( 1)-2),XA( 1),H,RRA,QQA,SA,AA,QA) 

END IF 

-CALCULATE CHANGE IN VARIABLES ACROSS SHOCK 

IF(SHOCK.EQ.322) THEN 
AB = AA/AR 

SAl s (Gl/G)i^OLOG( ( 2.0OO^^G^^(W^^it2)-G+l.DO0) / (G+l.DOO)) 
SA2 = Gli^DLOG( ( (G-l.DOO)^^(W**2) + 2.DOO) / ((G+l.DOO)^t 
C (Wit*2))) 

SB = SA4SA1>SA2 
UA » QA - VS 
UB = UA - ABJ^DQ 
QB = UB + VS 

ELSE 

CALL SKJUMP ( AA ,QA ,SA , AR ,DQ ,VS ,G ,G1 ,H , AB ,QB ,SB ) 

END IF 

•CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE 

IF(CNTACT.EQ.232) THEN 

QQ(I2(D) = QQCALC(QB,AB,SB) 

RR( 12(D) = RRCALC(QB,AB,SB) 

S(I2(D) = SB 
A(I2(D) = AB 
Q(I2(D) = QB 

ELSE 

QA s QB 
AA = AB 
SA = SB 

SB = S(I2(2)+1) 

CALL CSJUMP( AA,QA,SA,SB,G2,QB,AB) 

QQB = QQCALC(QB,AB,SB) 

RRB = RRCALC(QB,AB,SB) 

-INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES 

IF (12(2). EQ.N) THEN 

CALL BBDRY( RRB ,QQB ,SB ,RR( I2( 2 ) ) ,QQ( I2( 2 ) ) ,S( I2( 2 ) ) , 
A(I2(2)),Q(I2(2))) 

ELSE 

CALL INTERPC RRB ,RR( I2( 2 )♦! ) ,QQB,QQ( I2( 2 )♦! ) ,SB, 
S(I2(2)+1),XB(2),(XB(2)+H),RR(I2(2)), 
QQ(I2(2)),S(I2(2)),A(I2(2)),Q(I2(2))) 

END IF 
IF 



C 



C 

C 

END 

END IF 
GO TO ^0 
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BRANCH HERE IF LNCDE AND RNODE ARE CLOSE ENOUGH FOR SHOCK 

AND CONTACT SURFACE INTERACTION 

LEFT NODE: 

—SHOCK ON RIGHT, HEADED RIGHT, JUMPS NODE WITH CONTACT ON LEFT— 

—HEADED RIGHT, JUMPS NODE OR NOT > OR NO SHOCK WITH CONTACT ON — 

—LEFT, HEADED RIGHT, JUMPS NODE- 
RIGHT NODE: 

—SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE WITH NO CONTACT SURFACE- 

20 IF(( (LN0DE(2).EQ.221),AND,(LN0DE(3).EQ.321.0R,LN0DE(3).EQ.322)) 

C .0R.( LNODEt 2).EQ,100.AND.LNODE(3).EQ.32D) THEN 

IF( (RNODE(2).EQ.321).AND.(RNODE(3).EQ.100)) THEN 

DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE 

12(1) = RNODE(l) + 1 
IF(LNODE(2).EQ.322) THEN 
1212) = LNODEll) 

ELSE 

1212) = LNODEll) ♦ 1 

END IF 

CALL DELTAXl 12, SIGMA ,XB,XA,H) 

EXTRAPOLATE TO RIGHT FACE OF SHOCK 

IF 112(1). EG. N) THEN 

CALL BBDRYl RRl N ) ,QQ1 N ) ,S1 N ) ,RRB ,QQB ,SB ,AB,QB ) 

ELSE 

CALL EXTRAPl RRl 121 1 )), RRl I2( 1) + 1),QQ1 121 1) ),QQ1I2( !)♦!), 

C S(I2( 1) ),S1I2( 1)'M),XB1 1),H,RRB,QQB,SB,AB,QB) 

END IF 

CALCULATE VARIABLE CHANGE ACROSS SHOCK 

CALL SKJUMP(AB,QB,SB,AR,DQ,VS,G,G1,W,AA,QA,SA) 

DETERMINE CORRECTED VALUES AT NODEl 121 1 )-l ) AND ASSIGN THEM 

IF1LN0DE12).EG.100) THEN 

QQA 3 QQCALC1QA»AA,SA) 

RRA = RRCALC1QA,AA,SA) 

INTERPOLATE TO NODEl 121 1 )-l ) AND ASSIGN CORRECTED VALUES 

IF 112(1). EQ. 2) THEN 

CALL BBDRYl RRA ,QQA ,SA ,RR1 12( 1 )-l ) ,GQ1 12( 1 )-l ) , 

C S(I2( 1)-1),A(I21 1)-1),Q(I21 1)-D) 

ELSE 

CALL INTERPl RRA , RRl 121 1 )-2 ) ,QQA ,QQ1 121 1 )-2 ) ,SA , 

C S(I211)-2),XA11 ),1XA1 1)+H ),RR(I21 1)-1 ), 

C Q(31I21 1)-1),S(I21 1)-1) ,A1I21 1)-1), 

C (5112(1)-!)) 

END IF 

EXTRAPOLATE TO RIGHT FACE OF CONTACT SURFACE 

IF 11212). E(5.N) THEN 

CALL BBDRYl RRl N ) ,(5Q( N ) ,S( N ) ,RRB ,QQB ,SB , AB ,QB ) 

ELSE 

CALL EXTRAP1RR1I21 2)), RRl 12(2 )+l),C5Ql 12(2)), QQl 1212 )+l), 
C S(I212) ),S1I2(2) + 1),XB1 2),H,RRB,QQB,SB,AB,QB) 

END IF 

ELSE 

QQ(I2(1)-1) = Q(5CALC1QA,AA,SA) 

RR1I211)-1) = RRCALC1QA,AA,SA) 

S1I2(1)-1) = SA 
A1I211)-1) = AA , 

Q1I2(1)-1) = QA 
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IF(LNODE(3).EQ.522) THEN 

QQ(I2(2)) = QQ(I2(1)-1) 



RR(I2(21) 


= RR(I2(D- 


-1) 


S(I2(2)) 


= S( 


I2(l)-1) 


A(I2(21) 


» A( 


12(11- 


1) 


Q(I2(2)) 


= Q( 


12(1)- 


1) 


GO TO 40 









ELSE 

QB = QA 
AB s AA 
SB » SA 

END IF 

END IF 

SA = S(I2(2)-21 

CALL CSJUMP(AB>QB,SB,SA,62>QA>AA ) 

QQA - QQCALC(QA>AA>SA) 

RRA 3 RRCALC(QA>AA,SA1 

INTERPOLATE TO NODE( I2( 2 )-l ) AND ASSIGN CORRECTED VALUES 



C 



C 

C 

c 



IF (I2C2).EQ,2) THEN 

CALL BBDRY( RRA>QQA,SA,RR( I2( 2 )*1 ) >QQ(I2( 2 )-l 
S( I2( 2 )-l ) ,A( I2( 2 )-l) ,Q( I2( 2 )-l ) ) 

ELSE 

CALL INTERP( RRA>RR( 12(2 )-2 ) ,QQA>QQ( I2( 21-2 1 >SA > 
S(I2(2)-2),XA(2 ),(XA( 2UH ),RR(I2( 21-1), 

QQ(I2( 21-1), S( 12(2 )-l),A( 12(2 1-1), 

Q(I2(2)-in 

END IF 

END IF 



LEFT NODEj 

—SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE HITH NO CONTACT- 
RIGHT NODE: 

—NO SHOCK, WITH CONTACT ON RIGHT, HEADED LEFT, JUMPS NODE — 

ELSE IF((LNODE(2KEQ. 321). AND,(LN0DE(5).EQ. icon THEN 

IF( ( RNODE ( 2 j . £Q . 100 K AND . ( RNODE( 3 1 . EQ . 231 ) ) THEN 

DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE 

12(11 = LNODE(l) + 1 
12(2) = RNODE(l) 

CALL CELT AX( 12, SIGMA, XB,XA,H) 

CALCULATE JUMP THROUGH CONTACT SURFACE THEN SHOCK 

IF(I2C1).EQ.CI2(2)-11) THEN 
AA = A(I2(11) 

QA = Q( 12(111 
SA = S( 12(1 11 
SB = S( 12(2 1+1) 

CALL CSJUMP(AA,QA,SA,SB,G2,QB,AB) 

QQB 3 QQCALC(QB,AB,SB1 
RRB = RRCALC(QB,AB,SB) 

INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES 

IF (12(2). EQ.N) THEN 

CALL BBDRY(RRB,QQB,SB,RR(I2(2)),QQ(I2(2)),S(I2(2)), 
C A(I2(2)),Q(I2(2))1 

ELSE 

CALL INTERP( RRB,RR( I2( 2 )+l ) ,QQB,QQ( I2( 2 1+1 1 ,SB, 

C S(I2( 2) + l),XB(2 ),(XB(2)+H),RR(I2(21), 

C QQ(I2(2)),S(I2(21),A(I2(2)),Q(I2(2) )1 

END IF 

AB = A(I2(in 

QB - Q(I2(11) 

SB = S( 12(11) 
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CALL SKJUMP( AB,QB,SB,AR,DQ,VS,G,G1,H,AA,QA,SA) 
QQA « QQCALC(QA,AA,SA) 

RRA « RRCALC(QA>AA,SA) 

INTERPOLATE TO NOOE(I2m-l) AND ASSIGN CORRECTED VALUES 

IF a2(l),EQ.2) THEN 

CALL BB0RY( RRA ,QQA,SA,RR( I2( 1 )-l ) ,QQ( I2( 1 )-l) , 

C S(I2(1)-1),A(I2(1)-1),Q(I2(1)-1)) 

ELSE 

CALL INTERP( RRA, RR( 12(1 )-2), QQA, QQ(I2(1)-2),SA, 
C S(I2( 1)-2),XA( 1),(XA(1)+H),RR(I2( 

C QQ(I2(1)-1),S(I2(1)-1),A(I2(1)-1), 

C Q(I2(l)-in 

END IF 

ELSE 

RR(I2(1)-1) * RR(I2(l)-2) 

QQ(I2m-l) = QQ(I2(l)-2) 

S(I2(1)-1) = S(I2(l)-2) 

A(I2(1)-1) = A(I2(l)-2) 

Q(I2(1)-1) * Q(I2(l)-2) 

RR(I2(2)) = RR(I2(2) + 1) 

QQ(I2(2)) = QQ(I2(2)+1) 

S(I2(2)) = S(I2(2) + 1) 

A(I2(2n = A(I2(2) + 1) 

Q(I2(2)) = Q(I2(2) + 1) 

END IF 

END IF 



RIGHT NODE: 

—SHOCK ON RIGHT, HEADED LEFT, JUMPS NODE WITH NO CONTACT 

LEFT NODE: 

—-NO SHOCK, WITH CONTACT ON LEFT, HEADED RIGHT, JUMPS NODE 

ELSE IF( (RNODE( 2).EQ.231).AN0.(RNODE(3).EQ.100)) THEN 

IF( ( LNODE( 2 ) . EQ . 100 ) . AND . ( LNODE ( 3 ) . EQ . 321 ) ) THEN 

DETERMINE NODE TO RIGHT TO SHOCK AND CONTACT SURFACE 

12(1) = RNODE(l) 

12(2) = LNODE(l) + 1 

CALL 0ELTAX( 12, SIGMA, XB,XA,H) 

CALCULATE JUMP THROUGH CONTACT SURFACE THEN SHOCK 

IF(I2(2).£Q.(I2(1)-D) THEN 
AB = A(I2(2)) 

QB = Q(I2( 2) ) 

SB * S(I2( 2)) 

SA = S(I2( 2)-2) 

CALL CSJUMP( AB,QB,SB,SA,G2,QA,AA) 

QQA s QQCALC(QA,AA,SA) 

RRA s RRCALC(QA,AA,SA) 

INTERPOLATE TO NODE( I2( 2 )-l ) AND ASSIGN CORRECTED VALUES 

IF (I2(2),EQ.2) THEN 

CALL BBDRY( RRA,QQA,SA ,RR( I2( 2 )-l ) ,QQ( I2( 2 )-l ) , 
C S(I2( 2)-l),A(I2( 2)-l),Q(I2( 2)-D) 

ELSE 

CALL INTERP(RRA,RR(I2( 2)-2),QQA,QQ(I2( 2)-2),SA, 

C S(I2(2)-2),XA(2),(XA(2)+H),RR(I2(2)-1), 

C QQ( I2( 2 )-l ) ,S( I2( 2 )-l ) ,A( I2( 2 )-l ) , 

C Q(I2(2)-D) 

END IF 

AA = A(I2(2)) 

QA = Q(I2( 2) ) 

SA = S(I2(2)) 

CALL SKJUMP( AA,QA,SA,AR,DQ,VS,G,G1,W,AB,QB,SB) 
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QQ& B QQCALC(QB»AB>S&) 

RRB « RRCALC(QB,AB»SB) 

INTERPOLATE TO NODEC 12(D) AND ASSIGN CORRECTED VALUES 

IF (12(1). EQ.N) THEN 

CALL BBDRY( RRB>QQB,SB>RR( 12(D) >QQ( I2( 1) ) ,S( I2( D ) > 
C A(I2(1))>(3(I2(D)) 

ELSE 

CALL INTERP( RRB,RR( I2( 1 Ul ) ,QQB>(^( I2( 1 )4>1 ) ,SB, 

C S(I2( D + D,XB(D,(XB( 1)+H),RR(I2(D ), 

C (^(12(1) )tS(I2( 1))>A(I2( D)»Q( 12(D)) 

END IF 

ELSE 

RR(I2(2)-1) = RR(I2(2)-2) 

QQ(I2(2)-1) = QQ(I2(2)-2) 

S(I2(2)-D s S(l2(2)-2) 

A(I2(2)-D = A(I2(2)-2) 

Q(I2(2)-D = Q(I2(2)-2) 

RR(I2(D) = RR(I2(D+D 
QQ(I2(D) = QQ(I2(D + D 
S(I2(D) = S(I2(D + D 
A(I2(D) = A(I2(D+D 
Q(I2(D) = (3(I2(D + 1) 

END IF 

END IF 



RIGHT NODE: 

-—SHOCK ON LEFT, HEADED LEFT, JUMPS NODE HITH CONTACT ON RIGHT— 

HEADED LEFT, DOES NOT CROSS NODE OR NO SHOCK NITH CONTACT ON— 

RIGHT, HEADED LEFT, CROSSED NODE- 

LEFT NODE: 

SHOCK ON RIGHT, HEADED LEFT, JUMPS NODE WITH NO CONTACT— 

ELSE IF(( (RNODE(2).EQ.53D.AND,(RNO0E(3).EQ.231.OR.RNODE(5).EQ. 

C 23 2 ) ) . OR , ( RNODE ( 2 ) , E(3 . 100 , AND . RNOOE ( 3 ) . EQ . 231) ) THEN 

IF( ( LNODE( 2 ) .EQ. 231), AND. ( LNOOE( 3 ),EQ. 100 ) ) THEN 

DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE 

12(1) = LNOOE(l) 

IF(RN00E(3),EQ.232) THEN 

12(2) = RNODE(l) + 1 

ELSE 

12(2) = RNOOE(l) 

END IF 

CALL 0ELTAX( 12, SIGMA, XB,XA,H) 

EXTRAPOLATE TO LEFT FACE OF SHOCK 

IF (12(1). E(3, 2) THEN 

CALL BB0RY(RR( D,Q(3( D,S( D,RRA,QQA,SA,AA,QA) 

ELSE 

CALL EXTRAP(RR(I2( 1)-1),RR(I2(1)-2),QQ(I2(D-1),(3Q(I2( D-2), 
C S(I2( D-1),S(I2(D-2),XA(D,H,RRA,QQA,SA,AA,QA) 

END IF 

CALCULATE VARIABLE CHANGE ACROSS SHOCK 

CALL SKJUMP(AA,QA,SA,AR,0Q,VS,G,G1,W,AB,QB,SB) 

DETERMINE CORRECTED VALUES AT NOOE( 12(D) AND ASSIGN THEM 

IF(RN0DE(2).EQ.100) THEN 

QQB = QQCALC(QB,AB,SB) 

RRB = RRCALC(QB,AB,SB) 

INTERPOLATE TO NODE(I2(D) AND ASSIGN CORRECTED VALUES 

» 



183 



non ooooo ooo o o o ooo 



IF (12(1). EQ.N) THEN 

CALL BBDRY( RRB ,QQB ,SB,RR( 12(1)1 ,(JQ( I2( 1 ) ) ,S( I2( 1 ) ) , 
C A(I2(1))>Q(I2(1))) 

ELSE 

CALL INTERP(RRB,RR(I2(l)+l),QQB,OQ(I2(l)+l)fSBf 
C S(I2( l)4l)»XB(l)>(XB( l)4H)>RR(I2(l))f 

C QQ(I2(1)),S(I2(1) ),A( 12(D), Q( 12(1))) 

END IF 

EXTRAPOLATE TO LEFT FACE OF CONTACT SURFACE 



C 



IF (12(2). EQ. 2) THEN 

CALL BBDRY( RR( 1 ) ,QQ( 1 ) ,S( 1 ) ,RRA ,QQA ,SA ,AA,QA ) 

ELSE 

CALL EXTRAP( RR( I2( 2 )-l ) ,RR( I2( 2 )-2 ) ,QQ( I2( 2 )-l ) ,QQ( I2( 2 )-2 
S(I2( 2)-l),S(I2( 2)-2),XA( 2),H,RRA,QQA,SA,AA,QA) 

END IF 
ELSE 



QQ(I2(D) a QQCALC(QB,AB,SB) 
RR(I2(D) a RRCALC(QB,AB,SB) 
S(I2(D) a SB 
A(I2(D) a ab 
( 3(12(1)) a QB 
IF(RNODE(3).EQ.232) THEN 

RR(I2(1) + 1) a RR(I2(D) 
(3Q(I2(l)4l) a QQ(I2(D) 
S(I2(l)4l) = S(I2(in 
A(I2(l)4l) = A(I2(D) 
(3(12(1)41) s Q(I2(D) 
GO TO ^0 

ELSE 

(3A a QB 
AA a ab 
SA a SB 

END IF 

END IF 

SB a S(I2(2)4l) 



CALCULATE VARIABLE CHANGE ACROSS CONTACT SURFACE 

CALL CSJUMP( AA,QA,SA,SB,G2,QB,AB) 

(3QB a QQCALC(QB,AB,SB) 

RRB a RRCALC(QB,AB,SB) 

•INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES 



IF (12(2). EQ.N) THEN 

CALL BBDRY(RRB,QQB,SB,RR(I2(2) ),QQ(I2( 2)),S(I2(2))f 
C A(I2(2)),Q(I2(2))) 

ELSE 

CALL INTERP(RRB,RR(I2( 2)41), QQB,QQ(I2( 2)41), SB, 

C S(I2( 2)4l),XB( 2),(XB(2)4H),RR(I2( 2) ), 

C QQ(I2(2)),S(I2(2)),A(I2(2)),Q(I2( 2) )) 

END IF 

END IF 

END IF 
GO TO ^0 



BRANCH HERE IF THERE ARE ONE OR TWO NODES TO CORRECT, BUT THEY — - 

ARE SEPERATED BY MORE THAN 2H. CHECK FOR SHOCK/CONTACT SURFACE 

INTERACTION, AND CORRECT APPROPRIATELY — 



30 NODE a LNODE(l) 
SHOCK a LNODE(2) 
CNTACT a LN0DE(3) 



NO SHOCK, CONTACT SURFACE ON LEFT, HEADED RIGHT, JUMPS NODE 

35 IF( (SHOCK. EQ. 100). AND. (CNTACT. EQ. 321)) THEN 
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12(2) = NODE ♦ 1 

CALL 0ELTAX(I2>SI6t1A,XB>XA>H) 

CHECK FOR A SHOCK INTERACTING WITH CONTACT SURFACE AT NODE 

IF( ( SIGMA( 1 , 2 ) , LT. ( X2( 2 )4H ) ) . AND . 

C ( SIGMA( 1 , 2 ) . ST . ( X2 ( 2 ) ) ) )THEN 

AB = A(I2(2)) 

QB = (5(12(2)) 

SB = S(I2(2)) 

ELSE 

EXTRAPOLATE TO RIGHT FACE OF CONTACT SURFACE 

IF (12(2). HQ, N) THEN 

CALL BBDRY(RR(N)>QQ(N)>S(N),RRB»QQB»SB,AB,QB) 

ELSE 

CALL EXTRAP(RR(I2(2)),RR(I2(2)4l),QQ(I2(2)),QQ(I2(2)+l)> 
C S(I2( 2) ),S(I2(2 ) + l),XB( 2),H,RRB,QQB,SB,AB,QB) 

END IF 

END IF 

SA = S(I2(2)-2) 

CALCULATE VARIABLE CHANGE OVER CONTACT SURFACE 

CALL CSJUHP( AB,QB,SB>SA>G2>QA>AA) 

QQA s QQCALC(QA>AA,SA) 

RRA = RRCALC(QA>AA>SA) 

INTERPOLATE TO NODE( I2( 2 )-l ) AND ASSIGN CORRECTED VALUES 

IF (12(2). EQ, 2) THEN 

CALL BBDRY( RRA ,QQA ,SA ,RR( I2( 2 )-l ) ,QQ( I2( 2 )-l ) , 

C S(I2(2)-1),A(I2(2)-1),Q(I2( 2)-D) 

ELSE 

CALL INTERP( RRA»RR( I2( 2 )-2 ) ,QQA,QQ( I2( 2 )-2 ) >SA> 

C S(I2(2)-2),XA( 2),(XA( 2)+H),RR(I2( 2)-l), 

C QQ(I2( 2)-l),S(I2(2)-l),A( I2(2)-l), 

C Q(I2(2)-D) 

END IF 

— NO SHOCK, CONTACT SURFACE ON LEFT, HEADED LEFT, JUMPS NODE — 

ELSE IF( (SHOCK. EQ. 100 ). AND. (CNTACT.EQ. 231)) THEN 
12(2) = NODE 

CALL DELTAX( 12, SIGMA, XB,XA,H) 

CHECK FOR A SHOCK INTERACTING WITH CONTACT SURFACE AT NODE 

IF( ( SIGMA( 1 , 2 ) . GT , ( X2( 2 )-( 2 . DOO^^H ) ) ) . AND . ( SIGMA( 1 , 2 ) , LT . 

C (X2(2)-H))) THEN 

AA = A(I2(2)-1) 

QA = Q(I2(2)-1) 

SA = S(I2(2)-1) 

ELSE 

EXTRAPOLATE TO LEFT FACE OF CONTACT SURFACE 

IF (12(2). EQ. 2) THEN 

CALL BBDRY( RR( 1 ) ,QQ( 1 ) ,S( 1 ) ,RRA ,QQA ,SA, AA ,QA ) 

ELSE 

CALL EXTRAP( RR( I2( 2 )-l ) ,RR( I2( 2 )-2 ) ,QQ( I2( 2 )-l ) ,QQ( I2( 2 )-2 ) , 

C S( I2( 2 )-l ) ,S( I2( 2 )-2 ) ,XA( 2 ) ,H ,RRA ,QQA ,SA,AA ,QA ) 

END IF 
END IF 

SB = S(I2(2)4l) 

CALL CSJUMP( AA,QA,SA,SB,G2,QB,AB) 

QQB a QQCALC(QB,AB,SB) 

RRB = RRCALC(QB,AB,SB) 
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C INTERPOLATE TO NOOE(I2(2)) AND ASSIGN CORRECTED VALUES - 

C 

IF (I2t2),EQ.N) THEN 

CALL BBDRYC RRB ,QQB >SB ,RR( I2( 2 ) ) >QQ( I2( 2 ) ) >S( I2( 2 ) ) > 

C A(I2(2))>Q(I2(2))) 

ELSE 

CALL INTERPl RRB>RRt I2( Z)*l ) ,QQB,QQ( X2( Z)*l ) >SB> 

C S(I2t 2)+l),X8(2)>0<B(2)+H),RR(I2(2)), 

C QQ(I2(2) )>S(I2( 2n>A(I2(2n>Q(I2(2))) 

END IF 
C 

C SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE HITH NO CONTACT 

C 

ELSE IF( J SHOCK, EQ, 321), AND. tCNTACT.EQ, 100 n THEN 

12(1) = NODE + 1 

CALL DELTAX(I2, SIGMA, XB,XA,H) 

C 

C CHECK FOR A CONTACT SURFACE INTERACTING WITH SHOCK AT NODE 

C 

IF( (SIGMA( 2,2),LT,(X2( l)+H)),ANO. 

C (SIGMA(2,2),GT,(X2(1)))) THEN 

AB s A(I2(D) 

QB s Q(I2(1)) 

SB a S(I2(D) 

ELSE 

C 

C EXTRAPOLATE TO RIGHT FACE OF SHOCK 

C 

IF 112(1), EQ.N) THEN 

CALL BBDRY(RR(N),QQ(N),S(N),RRB,(^B,SB,AB,QB) 

ELSE 

CALL EXTRAP(RR(I2(1)),RR1I2(1)41),QQ(I2(1 )),QQ(I2( 1) + 1), 
C S(I2( 1) ),S( I2( 1 )4l) ,XB( 1),H,RR3,QQB,SB,A3,QB) 

END IF 

END IF 
C 

C CALCULATE VARIA3LE CHANGE ACROSS SHOCK 

C 

CALL SKJUMPl AB,QB,SB,AR,DQ,VS,G,G1,W,AA,QA,SA) 

QQA a QQCALC(QA,AA,SA ) 

RRA = RRCALC(QA,AA,SA) 

C 

C INTERPOLATE TO NODE! I2( 1 )-l ) AND ASSIGN CORRECTED VALUES 

C 

IF (12(1). EQ, 2) THEN 

CALL BBORYl RRA ,QQA ,SA ,RR( I2( 1 )-l ) ,QQ( I2( 1 )-l ) , 

C S(I2( 1 )-l),A(I2(l )-l),Q(I21 1)-1 )) 

ELSE 

CALL INTERPI RRA ,RR( I2( 1 )-2 ) ,QQA,QQ( I2( 1 )-2 ) ,SA, 

C S(I2( 1)-2),XA(1),(XA( 1)+H),RR(I2( 1)-1), 

C QQ(I2(1)-1),S(I2(1)-1),A(I2(1)-1), 

C Q(I2(1)-D) 

END IF 
C 

C — SHOCK ON RIGHT, HEADED LEFT, JUMPS NODE WITH NO CONTACT SURFACE— 

C 

ELSE IF( (SHOCK. EQ. 231). AND. (CNTACT.EQ, 100)) THEN 
12(1) = NODE 

CALL DELTAXl 12, SIGMA ,XB,XA,H) 

IF( (SIGMA(2,2).GT,(X2(1)-(2,D00J^H))),AND.(SIGMA( 2,2),LT. 

C (X2(l)-H))) THEN 

AA = A(I2(1)-1) 

QA a Q(I2(1)-1) 

SA a S(I2(1)-1) 

ELSE 

C 

C EXTRAPOLATE TO LEFT FACE OF CONTACT SURFACE 

C 
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IF (12(1). EQ. 2) THEN 

CALL EBDRY( RR( X ) fQQi X ) >S( 1 ) >RRA»QQA»SA,AA»QAi 

ELSE 

CALL EXTRAP(RR(I2(X)-X),RR(I2(X)-2),QQ(I2(X)-X),(5(5(I2(X)-2), 
C S(I2( X)-X)>S(I2(X)-2)>XA(X)>H9RRA>QQA»SA»AA>QA) 

END IF 
END IF 
C 

C — - CALCULATE VARIABLE CHANGE OVER SHOCK 

C 

IF(I2(X).EQ.N-X) THEN 
END IF 

CALL SKJUMP( AA,QA,SA,AR>DQ,VS,6,GX,H,AB>QB>SB) 

C 

C DETERMINE CORRECTED VALUES AT NODE(I2(X)) AND ASSIGN THEM 

C 

IF(I2(X).EQ.N-X) THEN 
END IF 

QQB s QQCALC(QB»AB»SB) 

RRB = RRCALC(QB,AB,SB) 

C 

c INTERPOLATE TO NODE(I2(X)) AND ASSIGN CORRECTED VALUES 

C 

IF (I2(X).EQ.N) THEN 

CALL BBDRY(RRB>QOB»SB,RR(I2(X) )>QQ(I2(X))»S(I2(X))» 

C A(I2(X)),(3(I2(X))) 

ELSE 

CALL INTERP( RRB,RR( I2( X )+X ) ,QQB ,QQ( 12( X )>X ) ,SB, 

C S(I2( X)+1),XB( X) ,(XB( X)+H) ,RR(I2(X ) ), 

C QQ( I2( 1) ) ,S( I2( 1) ) ,A( I2( X ) ) ,Q( I2( X ) ) ) 

END IF 

END IF 
C 

C CHECK IF SECOND NODE NEEDS CORRECTING IF SO LOOP BACK AROUND — 

C 



40 

C 

C 

C 

C 

C 

c 

c 

c 



c 

C CALCULATE DISTANCE NODE TO RIGHT OF DISCONTINUITY IS FROM LEFT — 

C BOUNDARY OF TUBE THEN DETERMINE DISTANCE FROM THIS NODE I2(*) 

C TO DISCONTI^>IUITY AND DISTANCE FROM NODE TO LEFT OF DISCONTINUITY- 

C TO THAT DISCONTINUITY 

C 

X2(X) = DBLE(I2(X)-X) * H 
XB(X) = (X2(X)-SIGMA(1.2)) 

XA( 1) = (H-XB( 1)) 

C 

C CALCULATE SAME VALUES FOR CONTACT SURFACE 

C 



IF(RNODE(X).EQ.O) GO TO 40 
NODE = RNODE(X) 

SHOCK = RN00E(2) 

CNTACT = RNOOE(3) 

RNODE(X) = 0 
GO TO 35 

RETURN 

END 

SUBROUTINE DELTAX( I2,SIGMA,XB,XA,H ) 

* LOCATE DISCONTINUITY WITHIN INTERVAL SUBROUTINE ^ 

^ # 

DIMENSION 12(4) ,SIGMA( 4,2) ,XB( 4 ) ,XA( 4 ) ,X2( 4 ) 

INTEGER 12 

DOUBLE PRECISION SIGMA ,XB,XA ,H ,X2 

X2(3)=O.DOO 

X2(4)=0.D00 
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C 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



X2(2) « DBLE(I2( 2)-l) H 
XB(2) s (X2(2)-SIGMA( 2,2)1 
XA(21 = (H-XB(21) 

RETURN 

END 

SUBROUTINE SKJUMP ( ADN ,QDN ,SDN , ARATIO , DELTAQ , VSHOCK ,G ,G1 , W , 
C AUP,QUP,SUP) 

w ww w W W WWWWW WWWW WW w w wwww w www w W WWW.WW WW W WWWW WW WW^W WW W W W W W WWW. 

SHOCK JUMP SUBROUTINE » 

VARIABLE DEFINITIONS 

*DN - VARIABLE DOWNSTREAM OF SHOCK 
ARATIO - SPEED OF SOUND RATIO 
DELTAQ - VELOCITY JUMP ACROSS SHOCK 
^\JP - VARIABLE UPSTREAM OF SHOCK 
VSHOCK - SHOCK VELOCITY 

DOUBLE PRECISION ADN, QDN, SDN, ARATIO, DELTAQ, VSHOCK, G,G1,W, 

C AUP ,QUP ,SUP ,SA1 ,SA2 ,UDN ,UUP 



•CALCULATE THE SPEED OF SOUND CHANGE THRU SHOCK 



AUP = ADN»ARATIO 



•CALCULATE THE ENTROPY CHANGE THRU SHOCK 



SAl = (Gl/G)i<DLOG((2.DOO^^G)f(W*i^2)-G+l.DOO)/(G+l.DOO)) 
SA2 = G1#DL0G( ( (G-l.D00)^t(W*^^2)^2)/( (G+l.D00)*(Wit^t2) )) 
SUP = SDN - SAl - SA2 



CALCULATE THE VELOCITY CHANGE ACROSS SHOCK 

UDN = QDN - VSHOCK 
UUP = UDN + ADN#DELTAQ 
QUP = UUP + VSHOCK 
RETURN 
END 

SUBROUTINE CSJUMP( ADN , QDN , SDN ,SUP ,G2 ,QUP ,AUP ) 

ik ^ 

^ CONTACT SURFACE JUMP SUBROUTINE » 



DOUBLE PRECISION QUP , QDN , AUP , ADN, SUP , SDN, G2 

CALCULATE THE SPEED OF SOUND CHANGE ACROSS THE CONTACT SURFACE- 
NOTE; THE VELOCITY IS CONTINUOUS ACROSS A C.S. 

QUP = QDN 

AUP = ADN DEXP( (SDN-SUP1/G2) 

RETURN 

END 

SUBROUTINE EXTRAP ( RRL ,RRR ,QQL ,QQR ,SL ,SR , DX , DH , RRIWT ,QQINT ,SIMT 
C AINT,QINT) 

» EXTRAPOLATION SUBROUTINE ^ 
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DOUBLE PRECISION RRL,RRR,QQL,OQR,SL,SR,DX,DH,RRINT,QQINT,SINT, 
C Air4T,QINT 



-REIMAN VARIABLES ARE EXTRAPOLATED, AND THE CORRESPONDING- 
-VALUES FOR VELOCITY AND SPEED OF SOUND ARE CALCULATED 



RRINT = RRL - (OX ( RRR-RRL VDH ) 

QQINT 2 QQL - iDX ^ ( QQR-QQL )/DH ) 

SINT = SL • (OX * (SR-SU/DH) 

AINT = (QQINT-RRINT)/(2,D00*SINT) 

QINT = (QQINT+RRINT)/( 2.D00) 

RETURN 

END 

SUBROUTINE INTERP( RRL ,RRR ,QQL ,QQR ,SL , SR ,DX ,DH , RRINT ,(3QINT ,SINT , 
C AINT, QINT) 

* INTERPOLATION SUBROUTINE * 

* * 

DOUBLE PRECISION RRL ,RRR ,QQL ,QQR ,SL , SR, OX, DH , RRINT , QQINT , SINT , 

C AINT, QINT 



-REIMAN VARIABLES ARE INTERPOLATED, AND THEIR CORRESPONDING- 
-VALUES OF VELOCITY AND SPEED OF SOUND ARE CALCULATED 



RRINT = RRL • (DX » ( RRL-RRR )/DH ) 

QQINT = QQL - ( DX * ( QQL-QQR )/DH ) 

SINT = SL • (DX * (SL-SR)ZDH) 

AINT = (QQINT-RRINT )/( 2.D00^^SINT) 

QINT = (QQINT+RRINT )/2.000 
RETURN 
END 

SUBROUTINE DBURST(N,H,QQ,RR,S,G,G1,G2,DELT,I2,X2,W,AR,DQ,VS, 
«LHPRES,SIGMA,A,Q) 

JH^**##**9^ It It If If ****»»*)^ **** * 

It It 

n DIAPHRAM BURSTING SUBROUTINE ^ 

It It 

itititititititititititititititititititititititititititititititititititititititititititititititit 

INTEGER N ,I , Y , 12 , L ,SHKDIR ,CSDIR , LHPRES 

DIMENSION X2(4),RR(N),QQ(N),S(N),I2(4),SIGMA(4,2),A(N),Q(N) 
DOUBLE PRECISION X2 ,X ,H , AB ,SA ,SB , AA ,QA ,QB ,QQA ,QQB , RRA ,RRB , SIGMA , 
C RR,QQ,S,TS,H,DQ,AR,PR,G,G1,G2,VS,DELT,CSRMN,A, 

C Q,MREIMN ,DREMN ,EREIMN ,WH ,SAX ,SA2 ,SAP ,SBP ,XB ,XA 

+++4++++ LOCATING THE NODE TO THE RIGHT +♦+♦♦♦+♦+++ 



DO 10 L=l,4 

SIGMA(L,1)=SIGMA(L,2) 

Y=0 

X=H 

1=2 

11 IF ( .NOT.(Y.EQ.O)) GOTO 10 

IF (SIGMA(L,1).LT.X) THEN 
X2( L)=X 
I2(L)=I 
Y=1 
END IF 
X=X+H 
1 = 1+1 
GOTO 11 
10 CONTINUE 
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CORRECT FOR OIAPHRAH BURSTING 



AT TIME ZERO DETERMINE CORRECT SHOCK DIRECTION 

SHKDIR = 3 IS A SHOCK HEADED LEFT, AND SHKDIR = E IS SHOCK 

TRAVELING RIGHT 

IF(LKPRES.EQ.3) THEN 
SHKDIR = 3 
X2(l) = X2(l) - H 
12 ( 1 ) = 12 ( 1 ) - 1 
X2(2) = X2(2) - H 
12 ( 2 ) = 12 ( 2 ) - 1 
GO TO 20 

ELSE 

SHKDIR = Z 

END IF 

20 RRA=RR(I2(1)-1) 

RRB=RR(I2(1) ) 

QQA=QQ(I2(1)-1) 

QGB=Qq(I2(l)) 

SA=S(I2(1)-1) 

SB=S(I2(1) ) 

21 AB=(QQB-RRB)/(2.000*SB) 

AA=( QQA-RRA )/( 2 . DOO^SA ) 

IF(SHKDIR.EQ.3) THEN 

MREIMN = (RRB-RRA)ZAA 
DREMN = DABSt MREIMN) 

ELSE 

MREIMN = (QQA-QQB)/AB 
DREMN = MREIMN 

END IF 

ITERATE FOR PROPER VALUE OF W USING THE QUADRATIC FIT OF THE 

REIMAN VARIABLE CHANGE HITH V CURVE. NOTE LEFT MOVING SHOCKS 

ARE USED IN THESE EQUATIONS SINCE RRB-RRA/AA = -( QQA-QQB/AB ) 

100 WW=( 3 . 0396408D01-( ( DREMN+2 . 7574D00 )/0 . 286337D00 ) ) 
W=5.513299000-DSQRT(HW) 

DQ=2.D00^(H»W-1.D00 )/( H*( G+1 . DOO ) ) 

AR=DSQRT( 2.0O0*(G-l.DOO))^(l.DOO + ( (G-1.000)»H»W/2.D00) )» 

C (G?^G2^H>H-1.D00 ) )/( (G + l.DOO )*H ) 

PR=( 2.D00*G/(G + 1.D00 ) )*W*H-( (G-l.DOO )/(G + l.D00 ) ) 

0R=( (G-1.000 )*W^H+2.D00 )/( (G + l.DOO )*H*H) 

EREIMN=DQ + ( AR-1 . DOO )»G2-( AR»G1/G )»DLOG( PR»( DR»»G ) ) 

IF ( DABS(EREIMN-DABS( MREIMN )).LT.0.1D-5) GO TO 110 
DREMN = (DABS(MREIMTJ) - EREIMN) + DREMN 
GOTO 100 

DETERMINE S BEHIND SHOCK AND THEN CALCULATE THE RIEMANN 

VARIABLE CHANGE ACROSS THE CONTACT SURFACE 

110 SAl = (G1/G)*DLOG((2.DOO*G»(W^^»2)-G + 1.DOO)/(G + 1.DOO)) 

SA2 = Gl»DLOG(((G-l.DOO)»(W»»2) + 2)/((G + l.DOO)»(H**2n) 
IF(SHKDIR.EQ.2) THEN 



SAP 


= 


SB - 


SAl - SA2 


SBP 


= 


SAP 




SBP 


s 


SA - 


SAl - SA2 


SAP 


= 


SBP 





END IF 

103 IF(SHKDIR.EQ.2) THEN 

CSRMN =( ( DEXP( ( SBP-SA )/G2 ) )»( SA )-( SBP ) )^^AR 

ELSE 

CSRMN =( ( DEXP( ( SAP-SB )/G2 ) )»( SB )-( SAP ) )»AR 

END IF 

IF(SHKDIR.EQ.3) THEN 
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MREIMN »((RRB-RRA)/AA)+CSRMN 
DREW4 = OABSJMREIMN) 

ELSE 

MREIMN \ QQA-QQB )/AB )-CSRMN 

OREMN = HREIMN 

END IF 

101 WH=( 5 . 0396408001-( ( OREMN+2 . 757^000 )/0. 286337000 ) ) 

W=5 * 51329^000-DSQRT ( WH ) 
0Q=2.000*(W)^Hrl.000)/(Wi^(G+l.D00)) 

AR=DSQRT( 2.D00i((G-l.D00 1 . D00+( ( G-1 . 000 )m^H/Z.DQQ ) 

C (G)^G2*WJ^W-1.000) )/( (G+1.D00)*W) 

PR=( 2 . D00^^G/( G+1 . 000 ) ( G-1 . 000 )/( G+1 . 000 ) ) 

0R=( (G-1. 000 ))^H*W+2.000 )/( ( G+1 . 000 ) 

EREIMN=0Q+( AR-1 . 000 ))^G2-( AR)^G1/G )^0L0Gi PR»( DR^^^^G ) ) 

IF (0ABS(EREIMN-DABS(MREIMN) ).LT. 0.10-5) GO TO 102 
OREMN = (OABS(MREIMN) - EREIMN) + OREMN 
GOTO 101 

102 SAl = (Gl/G)if0LOG((2.000*Gi((Wi^^^2)-G+l,000)/(G+1.000)) 

SA2 = Gl*0LOG(((G-1.000)*(WiH^2) + 2)/((G+1.000)^^(W^^*2))) 
IF(SHK0IR.EQ.2) THEN 

SAP = SB - SAl - SA2 

ELSE 

SBP = SA - SAl - SA2 

END IF 

IF (DABS(SAP-SBP).LT* 0.10-5) GO TO 105 
IF(SHKDIR.EQ.2) THEN 

SBP = (SAP-SBP) ♦ SBP 

ELSE 

SAP = (SBP-SAP) + SAP 

END IF 
GO TO 103 

CALCULATE SHOCK VELOCITY 

105 IF(SHKDIR.EQ.3) THEN 

DQ = ( -1.000 )^(DQ 

VS = ((RRA+QQA))^0,S000) - (W*AA) 

ELSE 

VS=( QQB+RRB )^^0 .S000+H*AB 

END IF 

IF(SHK0IR.EQ.2) THEN 

EXTRAPOLATE TO RIGHT FACE OF SHOCK 

CALL EXTRAP(RR(I2(1)),RR(I2( 1) + 1),(?Q( 12(1) ),QQ( 12(1 ) + l), 
C S(I2(1) ),S(I2( 1) + 1),H,H,RRB,QQB,SB>AB,QB ) 

CALCULATE CHANGE IN VARIABLES ACROSS SHOCK 

CALL SKJUMP(AB,(3B,SB,AR,DQ,VS,G,G1,W,AA,QA,SA) 

CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE 

QB = QA 
AB = AA 
SB = SA 

SA = S(I2(2)-2) 

CALL CSJUMP( AB>(3B>SB>SA>G2>QA»AA ) 

QQA 3 QA + AA)(SA 

RRA = QA - AA^^SA 

XA=0.00 

INTERPOLATE TO NOOE( I2( 2 )-l ) AND ASSIGN CORRECTED VALUES 

CALL INTERP( RRA ,RR( I2( 2 )-2 ) *QQA,QQ( I2( 2 )-2 ) ,SA, 

C S(I2( 2 )-2),XA,(XA+H),RR(I2( 2)-l), 

C QQ(I2(2)-1),S(I2(2)-1 ),A(I2(2)-1), 

C Q(I2(2)-1)>) 



ELSE 
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-EXTRAPOLATE TO LEFT FACE OF SHOCK 

CALL EXTRAP(RR(I2(l)-l),RR(I2(l)-2)>Qq(I2(l)-l)»QQCI2(l)-2)» 
C S(I2(l)-l),S(I2a)-2),H,H,RRA,QQA,SA,AA,QA) 

-CALCULATE CHANGE IN VARIABLES ACROSS SHOCK 

CALL SKJUMP(AA»QA>SA»AR>0Q»VS>G,G1>W>AB»QB»SB) 

-CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE 

QA = QB 
AA = AB 
SA = SB 

SB = S(I2(2)-H) 

CALL CSJUHP( AA>QA,SA,SB,G2>QB>AB) 

QQB = QB ♦ AB^SB 
RRB = QB - AB^SB 
XB=O.DO 



INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES 

CALL INTERP( RRB>RR( I2( 2 )♦! ) ,QQB,QQ( I2( 2 )-H),SB, 
C S(I2(2)-M )>XB,(XB>H)>RR(I2(2 n> 

C QQ(I2(2))»S(I2(2n>Aa2(2) )>Q(I2(2))) 

END IF 
RETURN 
END 



SUBROUTINE BONDRYC Q1 ,Q2 ,QB , A1 ,A2 ,QQ1 ,QQ2 ,RR1 ,RR2 ,S1 ,S2 ,H ,EE , 
C DELT ,BNDRY ,BDPRS ,BDPR ,BDDR ,BDTR , J , 

C NEWQQl ,NEWRR1 ,N£HS1 ,G ,G1 ,G2 ,HALT ,BDRY ,SK ) 

M BOUNDARY COUNDITION SUBROUTINE if 

»»»»»»»»»»»»»»» 



VARIABLE DEFINITIONS 

^^BD - VARIABLE AT PHANTOM NODE 

ifl - VARIABLE AT BOUNDARY NODE 

^^2 - VARIABLE AT FIRST NODE INSIDE BOUNDARY 

D2 - DENSITY RATIO AT FIRST NODE INSIDE BOUNDARY 

P2 - PRESSURE RATIO AT FIRST NODE INSIDE BOUNDARY 

QB - INITIAL VELOCITY AT AN OPEN BOUNDARY 

T2 - TEMPERATURE RATIO AT FIRST NODE INSIDE BOUNDARY 

DIMENSION AINT( 3 ) ,Z( 3 ) ,QPRIM( 3 ) , APRIMC 3 ) ,INTEG( 3 ) ,AAVG( 3 ) , 

C QQBD( 6 ) »RRBD( 6 ) ,SBD( 6 ) >QBD( 6 ) > ABD( 6 ) 

INTEGER BNDRY ,BDPRS > J ,HALT ,SK ,BDRY ,K , L ,M ,T 

DOUBLE PRECISION RRBD ,QQBD >SBD »QBD ,ABD ,BDPR ,BDDR >BDTR »DELQQH , 

C DELQQL ,DELRRH ,DELSH , DELSL ,DELQH , DELAH ,DELQL , 

C DELAL ,H ,EE , DELT ,QQINT ,RRINT ,SINT ,QPRIM, APRIM , 

C AINT,Z,SAVG,AAVG,Q1,Q2,A2,A1,RR1,QQ1,S1,QQ2,RR2, 

C S2,DLTAQQ,DLTARR,DLTAS,INTEG,QQSTEP,RRSTEP>D2, 

C SSTEP , NEWQQl ,NEWRR1 ,NEWS1 ,G ,G1 ,G2 ,DELRRL ,T2 , P2 , 

C AR »DQ , VS »W ,QB >NEWTR ,NEWDR ,NEWPR ,EE 1 ,EE2 >QRR >QQQ 

COMMON AR,OQ,VS,W 

- DETERMINE CURRENT PRESSURE AT NODE JUST PRIOR TO BOUNDARY 

T2 = ((QQ2-RR21»*2 V(^.D00»(S2»*2n 

D2 - ( ( l.D00/T2)»DEXP(Git( l,D00-G)*(S2-G2n )»^t(-Gl) 

P2 = T2 * D2 
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- DETERMINE IF LEFT OR RIGHT BOUNDARY 

IF(B0RY.EQ.2) THEN 
K a 1 
L a 2 
M a 5 

ELSE 

K a 6 
L a 5 
M a 4 

END IF 

IF( ( J.EQ.1).AN0,(BNDRY.EQ.0).AND.(B0PRS.EQ.0)) THEN 
QBD(K) a 0.000 
END IF 

QQBO(L) a QQi 
QQBD(M) a QQ2 
RRB0(U a RRl 
RRBD(M) a RR2 
SBD(L) a SI 
SBO(M) a S2 
QBO(L) a Ql 
QBD(M) a Q2 
AB0(U a A1 
ABD(M) a A2 

- DETERMINE IF BOUNDARY IS OPEN OR CLOSED, AND SET PROPER 

- VALUES AT PHANTOM NODE "LBD". BNDRY = 1( CLOSED ) ,0( OPEN ) 

S IF (BNDRY. EQ.l) THEN 

RRBOtK) a -QQBO(M) 

QQBO(K) a 0ABS(RRB0(M)) 

SBO(K) a SBO(M) 

QBO(K) a -QBO(M) 

ABO(K) a ABO(M) 

GO TO 10 

ELSE IF (BOPRS.EQ.l) THEN 
IF (J.EQ.l) THEN 

RRBD(K) a RRBD(L) 

QQBO(K) a QQBO(L) 

SBO(K) a SBO(L) 

ABO(K) a ABO(L) 

QBDiK) a QB 

END IF 
GO TO 10 

ELSE IF ((B0PR-P2).LE.0.1000) THEN 
IF(SK.E(3.3) GO TO 35 
S6D(K) a SBD(L) 

BOOR a ( (1.000/B0PR)^<(DEXPC(SB0(K)-G2)»(-(G/Gl)))n 
C »»(-(1.000/G)) 

BOTR a ((B00R)»i^(1.000/Gl))»0EXP(GiM1.000-G)»(SB0(K)-G2)) 
RRBO(K) a QBO( K )-( OSQRT( BDTR n^SBO( K I 
QQBO(K) a QB0(K) + (DSQRT(B0TR) l»SBD(K) 

ABO(K) a (QQBO(K) - RRBO(K)) / ( 2.D00»SBD( K ) ) 

RRBO(L) a RRBO(K) 

QQBO(L) a QQBO(K) 

SBO(L) a SBO(K) 

QBO(L) a QBD(K) 

ABO(L) a ABO(K) 

GO TO 10 

ELSE 

PRINT » ,'THE OPEN BOUNDARY HAS A PRESSURE HELD CONSTANT ’ 
PRINT ^ ,'THAT IS HIGHER THAN THE PRESSURE INSIDE THE TUBE* 
PRINT ^ , 'ADDITIONAL CODE IS NEEDED TO PRECEDE* 

, HALT a 1 
GO TO 20 

END IF 

-r CALCULATE THE JUMP TO THE NEXT TIME STEP USING CHARACTERISTICS— 
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C CHECK FOR SHOCK LEAVING TUBE OR HITHIN H OF BOUNDARY AND PICK — 

C CHOSE APPROPRIATE ALGORITHM TO ADVANCE BOUNDARY TO NEXT TIME ™ 

C 



10 IF(SK,HQ.2) THEN 

IF(B0RY.EQ.2) THEN 

CALL SKJUMP( ABD( L^l ) ,QBD( L>1 ) ,SBD( L4l ) ,AR,DQ»VS>G»G1 
C ABD(LnQBO(L)>SBD(D) 

QQBD(L) = QBD(U 4 ABD( L))^SBD( L ) 

RRBO(L) = QBO(U - ABD( L )^^SBO( L ) 

GO TO 35 

ELSE IF(BDRY.EQ.3) THEN 

CALL SKJUMPl ABO( L-1 ) ,QBD( L-1 1 ,SBO( L-1 ) ,AR ,DQ ,VS ,G,G1 ,W, 
C ABD(L)»QBD(L)>SBD(Ln 

QQBD(L1 = QBD(L) 4 ABD( U»^SBO( L ) 

RRBD(L) = QBD(U - ABD( L )^SBD( L 1 
GO TO 35 

END IF 

END IF 

OELQQL = QQBO(L) - QQBO(L-l) 

DELQQH = QQBD(L4l) - QQBO(L) 

OELRRH = RRBD(L4l) - RRBO(U 
DELRRL = RRBD(U - RRBD(L-l) 

DELSH = SBD(L4l) - SBD(U 
DELSL = SBD(L) - SBD(L-l) 

DELQH = QBD(L4l) - QBO(L) 

DELQL = QBO(L) - QBO(L-l) 

DELAH = ABD(Lil) - ABO( L ) 

DELAL = ABD(L1 - ABD( L-1) 

IF( (BNORY,EQ.O).AND.CSK.EQ.in THEN 
IF (BDRY.EQ.2) THEN 

CALL C0ND3( QBD( L ) ,QBO( L+l ) , ABD( L ) , ABD( L+1 ) ,RRBD( L ) ,QQBO( L ) , 

C SBO( L ), DELQQH , OELRRH , DELSH , DELQH , DELAH ,OELT ,H ,EE , 

C QQINT,RRINT,SINT,QPRIM,APRIM,AINT ) 

ELSE 

CALL CON02( QBD( L ) ,QBO( L-1 ) ,ABD( L ) ,ABO( L-1 ) ,RRBO( L ) ,QQBD( L ) , 

C SBO( L),DELQQL, DELRRL, DELSL, OELQL, DELAL, CELT, H, EE, 

C QQINT ,RRINT ,SINT ,QPRIM ,APRIM,AINT ) 

END IF 

ELSE IF((BNDRY,EQ.1).AN0.(SK.EQ.D) THEN 
QQSTEP = 0.000 
RRSTEP = O.DOO 
SSTEP = O.DOO 
GO TO 30 
C 

C USE CONDITION 1 ALGORITHM TO CALCULATE RRINT , QQINT , SINT 

C 

ELSE 

CALL CONOK QBD( L ) ,QBO( L + l ) , ABO( L ) ,ABD( L+1 ) ,RRBD( L ) ,QQBO( L ) , 



C SBO( L ) ,DELQQL, OELRRH , DELSH , DELSL , 

C DELQH, DELAH, DELQL, DELAL, H, EE, DELT, QQINT, RRINT, SINT, 

C QPRIM,APRIM,AINT,QBD(L-1 ),ABD( L-1 ) ) 

END IF 
C 

c CALCULATE OLTA QQ, DLTA RR i DLTA S 

C 



DLTAQQ=QQINT-QQBO( L ) 

DLTARR=RRINT-RRBD( L ) 

OLTAS=SINT-SBO(U 

C 

c CALCULATE Z(K)*S 

C 

C 

AAVG( 1 )=( AINT( 1 )+ABO( L ) )/2,0000 

AAVG( 3 ) = ( AINT( 3 )4ABD( L ) )/2 . 0000 

AAVG(2)=0.0D00 

SAVG = (SINT4SB0( U)/2.000 

Z( 1 )=-( 1 . 0000/G2 )^AAVG( 1 )^( SAVG-G2 )^( QPRIM( 1 )-G2^APRIM( 1 ) ) 
Z( 3 ) = ( 1 . ODOO/G2 ))(AAVG( 3 )^( SAVG-G2 QPRIM( 3 )4G2^APRIM( 3 ) ) 
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Z(2)=0.0000 

INTEGRATE THE Z(K)*S 

INTEG( 2) =0.0000 

INTEG(1)=Z(1)5^0ELT 

INTEG(3)=ZI3il^^0ELT 

SOLVE THE EQUATION 

QQSTEP=0UTAQQ+INTEG( 1 ) 

RRSTEP=DLTARR+INTEG( 3 ) 

SSTEP=DLTAS+INTEGI 2 ) 

STORE THE SOLUTION 

30 NEWQQ1=QQBD( L )+QQSTEP 

NEWRR1=RRBD( L ) +RRSTEP 
NEWS1=SBD( U+SSTEP 

UPOATE PHANTOM NODES FOR OPEN BOUNDARY CONDITION 



C 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



35 IF ((BNDRY.EQ.0).AND.(BDPRS.EQ.1))THEN 
RRBD(K) = RRBD(U 
QQBO(K) = QQBD(L) 

SBD(K) = SBD(L) 

ABD(K) = ABD(L) 

QBD(K) = QBD(L) 

END IF 

IF ( (BNDRY.EQ.O).AND.(BDPRS.EQ.On THEN 
IF (SK.EQ.2) THEN 

QBD(K) = (QQBD(U‘»RRBD(L)) / 2.000 
ELSE 

NEWTR= ((NEWQQl - NEWRRl )*J^2 )/( 4. D00J(( ( NEWSl )*^^2 ) ) 
NEWDR= ( ( 1 . DOO/NEWTR )»DEXP( G^^ ( 1 . DOO-G m NEWS1-G2 ) ) ) 

C 

NEWPR= NEWTRifNEWOR 
EEl = DABS(BDPR-NEWPR) 

EE2 = 0ABS(SBD(L)-NEWS1) 

IF (CEEl. GT. 0.10-5). OR. (EE 2 .GT. 0 . 1 D-Sn THEN 
QRR= NEWRRl + ( DSQRT( BDTR ) )*SBD( L ) 

QQQ= NEWQQl - ( DSQRT( BDTR ) )*SBD( L ) 

QBD(K) = (QRR+QQQ)/2.DDD 
GO TO 5 

END IF 
END IF 

END IF 

20 QQl = QQBO(L) 

RRl = RRBDtL) 

SI = SBD(L) 

Q1 = QBD(L) 

A1 = ABD(L) 

RETURN 

END 

SUBROUTINE SRFLCTC QQA ,RRA ,SA , SIGMA ,VS ,DELT , LWPRES ,RRB ,QQB ,SB ,QB , 
C AB,G,G1,G2) 

M * 

^ SHOCK REFLECTION AT SOLID BOUNDARY SUBROUTINE * 



VARIABLE DEFINITIONS 

OELTEX - THE EXCESS TIME IN A TIME STEP WHEN THE SHOCK IS 
EXACTLY AT THE SOLID WALL 
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OELTWL - THE TIME FOR THE SHOCK TO REACH THE WALL 

DIMENSION SI6MA(<f»2) 
irrrEGER lwpres 

DOUBLE PRECISION QA>QB»AA,AB>SA,SB»QQA>QQB,RRA,RRB »SI6HA, 

C DQ,W,AR,PR,DR,EREIMN,SA1,SA2,VS,DELT,DELTWL, 

C DELTEX,G>G1>G2 

CALCULATE THE TIME FOR SHOCK TO REACH WALL, AND EXCESS TIME — 

IN THIS TIME STEP 

IF (LWPRES. EQ. 2) THEN 
DELTWL = (1.000-SIGMA(l,in/VS 
DELTEX = DELT - DELTWL 

CALCULATE VELOCITY GRADIENT OVER SHOCK AS IT REFLECTS 

QA = (QQA>RRA)/2.D00 
QB = O.DOO 

AA = (QQA-RRA)/(2.D00*SA) 

DQ = (QB-QA)/AA 

CALCULATE W 

W = DSQRTt ( (DQ*^^2)»0.36D00) + 1.D00) - (DQ^^0.6D00) 

FOR SHOCK AT LEFT BOUNDRY DO THE SAME 

ELSE 

QQB=QQA 

RRB=RRA 

SB=SA 

DELTWL = (DABS(SIGMA(l,in-O.DOO)/DABS(VS) 

DELTEX = DELT - DELTWL 

CALCULATE VELOCITY GRADIEhfT OVER SHOCK AS IT REFLECTS 

QB = (QQB+RRB)/2.D00 
QA = O.DOO 

AB = (QQB-RRB)/(2.D00^^SB1 
DQ = (QA-QBl/AB 

CALCULATE EXACT REIMAN VARIABLE JUMP FROM DQ 

W = DSQRT( ( (DQ**2)^^0.36D00) + 1.D00) + (DQ*0.6D00) 

END IF 

CALCULATE AR,PR,DR OVER SHOCK 

AR=DSQRT( 2.D00^^(G-1.D00)^^( 1. DOO+( ( G-l.DOO )^^W*W/2. DOO ) )» 

C (G^^G2*W)^W-l.D00n/( (G+l.DOO )*W ) 

PR=( 2, DOO*G/( G+l.DOO ) )*W*W-(( G-1 . DOO )/( G+l.DOO)) 
DR=((G-1.D00)^^W*W+2.D00 )/( (G+l.DOO )^^W^W ) 

SAl = (Gl/G)?^DLOG( ( 2.D00J^G»(W*^^2)-G+1.000)/(G+1.D00) ) 

SA2 = GH^DLOG(((G-1.0O0)^(W^^*2) + 2)/((G+l.DOO)*(W^^^^2))) 

CALCULATE ENTROPY AND SPEED OF SOUND BEHIND SHOCK 

IF (LWPRES. EQ. 2) THEN 
SB = SA - SAl - SA2 
AB = AA)^AR 

CORRECT REIMAN VARIABLES AT BOUNDARY 

RRB = QB - ABJ^SB 
QQB = QB + AB>*SB 

CALCULATE NEW SHOCK SPEED 
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vs « l(RRA4QQA)«0.5D00) - (W«AA) 

CALCULATE POSITION OF REFLECTED SHOCK AT END OF THIS TIME STEP- 

SIGMA(1,2) = VS*DELTEX ♦ l.DOO 
LWPRES * 3 

SHOCK REFLECTING AT LEFT BOUNDARY 

ELSE 

SA » SB - SAl - SA2 
AA 3 AB^AR 

CORRECT REIMAN VARIABLES AT BOUNDARY 

RRA 3 QA - AA^SA 
QQA 3 QA AA^SA 

CALCULATE NEW SHOCK SPEED 

VS 3 ( (RRB+QQB)*0.5D00) ♦ ( W*AB ) 

CALCULATE POSITION OF REFLECTED SHOCK AT END OF THIS TIME STEP- 

SI6MAC1,2) 3 VSmDELTEX 
LWPRES = 2 
RRB3RRA 
QQB^QQA 
SB=SA 
QBsQA 
AB=AA 
END IF 
RETURN 
END 

SUBROCTTINE BBDRY( RRI >QQI >SI >RR2 >QQ2 >S2 > A2 >Q2 ) 

* * 

» NODES NEAR BOUNDARY SUBROUTINE * 

DOUBLE PRECISION RRI ,QQ1 ,S1,RR2,QQ2 ,S2 ,A2 ,Q2 
RR2 = RRI 
QQ2 3 QQl 
S2 3 SI 

A2 3 (QQ2 - RR2)/( 2.DOO*S2) 

02 3 (QQ2 > RR2)/2.000 
RETURN 
END 



SUBROUTINE B0RDER( JSTOP 1 

* # 
* PLOTTING AREA SETUP ROUTINE * 

INTEGER I 

REAL X0RGC^)/1.75,^.65,1.75,^.65/ 

REAL YORG ( ^ )/2 . 75 , 2 • 75 , 5 . 80 , 5 . 80/ 

REAL YMAX(4)/5.5,5.5,1.5,6.20/ 

REAL YMIN(4)/0.5,0.5,-0.5,4.90/ 

DO 70 1=1,4 

CALL PHYSOR(XORG(I),YORG(in 
CALL NOBROR 
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CALL AREA2D(2.4,2.4) 

CALL FRAME 

CALL GRAF( 0 . , * SCALE * ,1 . 0 ,YMIN( 1 ) , 'SCALE ' ,YMAX( I ) ) 

CALL ENOGR(O) 

70 CONTINUE 

CALL PHYSOR( 1.25,2. 25) 

CALL NOBROR 

CALL AREA20(6.00,6.50) 

CALL HEAOIN( 'SHOCK TUBE RESULTS$' ,100,1.7,4) 

CALL HEAOIN( 'FIRST ORDER N =101$ * ,100 ,1. 2 ,4 ) 

CALL HEA0IN( 'DENSITY RATIO = 5.0 TEMP RATIO = 1$ ' ,100 ,1 . 0 ,4 ) 
CALL HEADIN( 'PRESSURE RATIO = 5.0$' ,100,1.0,4 ) 

CALL GRAF(0. , 'SCALE' ,1. ,0. , 'SCALE' , JSTOP ) 

CALL ENOGR(O) 

RETURN 

END 



SUBROUTINE PLOT( J , JSTOP ,N ,QQ ,RR ,S ,H ,XARRAY , 
ffPARRAY ,DARRAY ,QARRAY ,SARRAY ,G ,G1 ,G2 ) 

it it 

it GRAPHICAL PLOTTING ROUTINE it 

it it 

ititititititititititititititititititititititititititititititititititititititititititititititit 

INTEGER I ,N , J , JSTOP ,KNT( 4 )/l ,4 , 6 , 9/ 

DIMENSION QQ( N ) ,RR( N ) ,S( N ) ,XARRAY( N ) ,QARRAY( N ) ,PARRAYt N ) , 

C DARRAY( N ) ,SARRAY( N ) 

DOUBLE PRECISION QQ ,RR ,S ,H ,G ,G1 ,G2 
REAL XARRAY,QARRAY,PARRAY,DARRAY 
REAL SARRAY,TEMP,HR,G1R,GR,G2R 
REAL XORG( 4 )/l . 75 ,4 . 65 , 1 . 75 ,4 . 65/ 

REAL YORG( 4 )/2 . 75 , 2 . 75 ,5 . 30 ,5 . 30/ 

REAL YMAX(4)/5.S,5.5,1.0,6.20/ 

REAL YMIN(4)/0.5,0.5,-1.0,4.90/ 

CHARACTER)t4 lYNAM 
DIMENSION lYNAMdS) 

DATA lYNAM/ ' PRES ' , ' SURE ' , ' $ ' , ' DENS ' , 

ff'ITY$','VELO','CITY','$ ' , 'MODI ' , ' FIED ' , ' ENT ', 'ROPY ',' $ '/ 

CONVERT DOUBLE PRECISION TO SINGLE PRECISION 

GR=SNGL(G) 

G1R=SNGL(G1) 

G2R=SNGL(G2) 

HR=SNGL(H) 

DO 80 1=1, N 

QARRAY( I ) = ( SNGL( QQ( I)+RR( I) )/2. 0 ) 

TEMP=SNGL( QQ( I )-RR( I ) )J^SNGL( QQ( I )-RR( I ))/( 4 . OJ^SNGL( S( I )^S( I) ) ) 
OARRAY( I ) = ( ( 1 . 0/TEMP )itEXP i GRit{ 1 . 0-GR )it{ SNGL( S( I ) )-G2R ) ) )ititi -GIR ) 
SARRAY(I)=SNGL(Sa) ) 

PARRAYl I )=TEMP*DARRAY( I ) 

80 CONTINUE 

DO 83 1=1,4 

CALL PHYSOR(XORG(I),YORG(D) 

CALL AREA2D( 2.4,2.4) 

CALL XNAME( 'X' ,1) 

CALL YNAME(IYNAM(KNT(I)),100) 

CALL GRAF( 0 . , ' SCALE ' , 1 . 0 ,YMIN( I ) , ' SCALE ' ,YMAX( I ) ) 

IF CI.EQ.l) CALL CURVE! XARRAY ,P ARRAY ,N ,0 ) 

IF (I.EQ.2) CALL CURVE! XARRAY ,DARRAY ,N ,0 ) 

IF (I.EQ.3) CALL CURVE! XARRAY ,QARRAY ,N,0 ) 

IF !I.EQ.4) CALL CURVE ! XARRAY, SARRAY, N, 0 ) 

CALL ENDGR!0) 

33 CONTINUE 
RETURN 
END 
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SUBRCXJTINE EXACT(N*XINir,T,VHEAO,VTAIL,VCDE,VSE>DLCO,OLSH,QQ,RR,S> 
nn >XARRAY >DARRAY yG »Q1 >62 >CRI ) 



i* EXACT SOLUTION COMPARISON ROUTINE » 

» » 



INTEGER N 

DIMENSION QQ( N ) ,RR( N ) ,S( N ) ,XEXACT ( 6 ) ,YEXACT( 6 ) ,XARRAY ( N ) , 

C DARRAY(N) 

DOUBLE PRECISION QQ, RR, S,H, 6,61 ,62,7, DRI 
REAL XEXACT ,YEXACT ,XINIT , VHEAD , VTAIL ,VCDE 
REAL VSE ,DLCD ,DLSH ,TEMP ,XARRAY >DARRAY 

XEXACTC 11=0.0 

XEXACT ( 2 1=XINIT+T»VHEAD 

XEXACT ( 3 )=XINIT+Ti^VTAIL 

XEXACK 4 )=XINIT+TJ^VCDE 

XEXACT ( 5 )=XINIT+T^tVSE 

XEXACT(6)=1.0 

YEXACT(1)=SNGL(DRI) 

YEXACT(2)=YEXACT(1) 

YEXACT(3)=DLCD 

YEXACT(4)=YEXACT(3) 

YEXACT(5)=DLSH 

YEXACT(6)=1.0 

DO 90 1=1, N 
QQ(I)=SNGL(OQ(in 
RR(I )=SNGL(RR(I n 
S(Il=SNGL(Sri n 

TEMP=( QQ( I )-RR( 1 1 )»( QQ( 1 1-RR( 111/(4. 0J^S( I )»S( I 1 1 
DARRAYC 1 1 = ( ( 1 . 0/TEMP l^^EXPC GM 1 . 0-6 1^( S( 1 1-62 1 1 -61 1 

90 CONTINUE 

CALL PHYSOR( 2.65,2.251 
CALL NOBRDR 
CALL AREA2D(6.75,4.01 
CALL XNAME( *X* ,11 
CALL YNAME( 'DENSITY* ,71 

CALL HEADIN( 'DENSITY OISTRIBUTION$ ' ,100 ,1.3 ,4 1 

CALL HEADIN( 'FIRST ORDER N =101$ ' ,100 ,1. 0 ,4 1 

CALL HE ADIN( 'DENSITY RATIO = 5 TEMP RATIO = 1$ * ,100 , 1.0 ,4 1 

CALL HEADIN( 'PRESSURE RATIO = 5 $ ' , 100 ,1 . 0,4 1 

CALL LINES( 'EULER-1 SOLUTION$ * ,IPKRAY ,2 1 

CALL LINES( 'EXACT SOLUTION$ * ,IPKRAY ,1 1 

CALL FRAME 

CALL GRAF( 0.0, 'SCALE* ,1.0, 0.0, 'SCALE' ,6.01 
CALL MARKER(01 

CALL CURVE(XEXACT,YEXACT,6,-11 
CALL LEGLIN 

CALL CURVE (XARRAY,OARRAY,N, 01 
CALL LEGEND(IPKRAY,2,4.25,2.75) 

CALL ENDGRCOl 

RETURN 

END 
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